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Abstract 


The  use  of  composite  material  systems  in  structural  design  has  yielded 
significant  improvements  in  material  efficiency  by  minimizing  weight  while 
meeting  static  and  dynamic  strength  requirements.  The  complexity  of  composite 
materials,  however,  has  presented  a  spectrum  of  additional  design 
considerations  in  the  areas  of  fabrication,  strength  tailoring,  failure  mechanisms, 
and  damage  tolerance.  Laminated  composites  exhibit  a  variety  of  damage  and 
failure  modes  which  include  matrix  cracking,  fiber  breakage,  delamination 
propagation,  and  instability.  Of  prime  concern  is  the  predication  of  residual 
strength  in  composite  structures  subjected  to  damage-causing  impact  threats  in 
the  intended  service  environment.  The  present  effort  aims  to  account  for  many 
of  the  salient  damage  mechanisms  utilizing  a  numerical,  finite  element-based 
approach.  Such  an  approach  permits  a  broad  range  of  failure  modes  to  be 
accounted  for  while  allowing  the  modeling  of  complex  geometries,  support 
conditions,  and  applied  loading.  This  report  details  the  theoretical  and 
algorithmic  basis  of  a  developed  computer  program  denoted  RESTRAN 
(REsidual  STRength  ANalysis),  incorporating  both  material  and  structural  failure 
modes  in  the  prediction  of  residual  strength  in  composite  structures  containing 
internal  damage. 


Table  of  Contents 


Page 


List  of  Figures  .  vii 

List  of  Tables  . ix 

1.  Introduction  .  1 

2.  Finite  Element  Formulation  .  2 

2.1  Coordinate  Systems  and  Nomenclature . 3 

2.1.1  Global  and  Local  Coordinates  .  3 

2.1.2  Coordinate  Transformations  .  6 

2.2  Elastic  Stiffness  Formulation  .  7 

2.3  Modeling  Material  Nonlinearity  .  11 

2.4  Modeling  Geometric  Nonlinearity  .  14 

3.  Laminated  Material  Representation  .  18 

3.1  Effective  Material  Moduli  .  18 

3.2  Ply-Level  Stress  and  Strain  Recovery  .  21 

3.3  Nonlinear  Elastic  Moduli  .  22 

4.  Material  Failure  Modes  .  24 

4.1  Material  Failure  Criteria .  24 

4.1.1  Maximum  Stress  Criterion  .  27 

4.1.2  Maximum  Strain  Criterion  .  27 

4.1.3  Beltrami  Criterion  .  28 

4.1.4  Von  Mises  Criterion  .  28 

4.1.5  Hoffman  Criterion . 29 

4.1.6  Hill  Criterion  .  29 

4.1.7  Tsai-Wu  Criterion . 30 

iii 


4.1.8  Christensen  Criterion  . 

4.1.9  Feng  Criterion  .  31 

4.1.10  Modified  Hashin  Criterion . 32 

4.2  Progressive  Material  Failure  Prediction  .  33 

4.3  Material  Damage  Models  . 33 

5.  Structural  Failure  Modes  . 35 

5.1  Multiple  Delaminations  .  37 

5.2  Algorithmic  Assessment  of  Local  Buckling  Failure . 39 

5.3  Distributed  Delaminations  and  Simultaneous  Local  Buckling  .  46 

5.4  Post-Buckled  Failure  Modeling . 47 

6.  Solution  Algorithms . 49 

6.1  Prepass . . 49 

6.2  Static  and  Buckling  Analysis . 50 

6.3  Residual  Strength  Analysis  Algorithm  .  53 

7.  Ultimate  Failure  Prediction  .  55 

8.  Computer  Implementation  of  RESTRAN .  55 

9.  Numerical  Studies . 57 

10.  Practical  Analysis  of  Residual  Strength  .  55 

10.1  Damage  Characterization  . . .  56 

10.2  Demonstration  Problem  .  57 

10.2.1  Input  Data  File  for  a  Residual  Strength  Problem  .  68 

10.2.2  Output  Data  File  for  a  Residual  Strength  Problem  .  77 

10.2.3  Graphical  Output  Using  MATHEMATICA  . 84 

11.  Conclusion . 38 


IV 


12.  References  .  89 

Distribution  List . 97 

Report  Documentation  Page  .  115 


v 


INTENTIONALLY  LEFT  BLANK. 


List  of  Figures 


Figure  Page 

I.  Global  and  local  coordinate  systems . 3 

2(a).  Laminate  layup  description .  4 

2(b).  Principal  ply  directions .  4 

2(c).  Ply  orientation  angle  convention .  5 

3.  Stress  and  strain  tensor  notation .  6 

4.  Local  (x’,y’,z’)  and  (x,y,  z)  Coordinate  Systems .  6 

5.  Hexahedral  element  configuration . 9 

6.  Partial  failure  in  nonlinear  elastic  material .  12 

7.  Iterative  solution  for  nonlinear  materials .  13 

8.  Bifurcation  of  equilibrium  states .  15 

9.  Nonlinear  elastic  stress-strain  relation .  22 

10.  Secant  moduli  depiction .  23 

II.  Global,  local  and  mixed  buckling  modes  in  thick  laminates .  36 

12.  Laminate  with  multiple  embedded  delaminations .  37 

13.  Example  of  a  nonphysical  sublaminate  buckling  mode .  37 

14.  Condensation  of  multiple  delaminations .  38 

15.  Sequential  buckling  failure .  39 

16.  Local  compatibility  violation  in  layer  buckling .  39 

17.  Designation  of  element  faces .  40 

18.  Generation  of  coincident  nodes . 40 

19.  Opening  and  closing  buckling  modes .  41 

20.  Nodes  used  to  calculate  normal  modal  strain .  42 

vii 


21.  Definition  of  positive  normal  to  delamination  plane .  43 

22.  Opening  mode  in  delamination  buckling . 44 

23.  Buckling  failure  in  a  centrally  located  delamination .  45 

24.  Failure  assessment  of  buckled  sublaminate .  45 

25.  Modeling  traversing  delaminations .  46 

26.  Arbitrarily  oriented  delaminations .  47 

27.  Post-buckling  load  carrying  capability .  48 

28.  RESTRAN  solution  control  options .  49 

29.  Prepass  flow  chart .  50 

30.  Flow  chart  of  static  analysis . 51 

31.  Flow  chart  of  buckling  analysis .  52 

32.  Flow  chart  of  residual  strength  solution  algorithm . 54 

33.  Models  for  Euler  column  buckling  load  determination .  58 

34.  Geometry  and  loading  of  a  flat  plate .  59 

35.  Laminated  plate  with  elliptical  delamination .  52 

36(a).  NDIV  =  4.  No  constraints  applied.  A  =  2.0525F75 .  62 

36(b).  NDIV  =  4.  Compatibilty  constraints  enforced.  A  =  3.1637.E5 . 63 

36(c).  NDIV  =  8.  No  constraints  applied.  A  =  1.4667J55 .  63 

36(d).  NDIV  =  8.  Compatibility  constraints  enforced.  A  =  2.1523E5 .  63 

36(e).  NDIV  =  12.  No  constraints  enforced.  A  =  1.3864U5 .  64 

36(f).  NDIV  =  12.  Compatiblity  constraints  enforced.  A  =  1.9803E5 . 64 

36(g).  NDIV  =  16.  No  constraints  applied.  A  =  1.3956E5 .  64 

36(h).  NDIV  =  16.  Compatibility  constraints  applied.  A  =  2.0527E5 . 65 

37.  Geometry  and  loading  of  an  elastically  supported  composite  plate . 67 

viii 


List  of  Tables 

Table  Page 

1.  Maximum  stress  criterion  failure  modes .  27 

2.  Maximum  strain  criterion  failure  modes .  28 

3.  RESTRAN  solution  options .  53 

4.  Tests  for  catastrophic  failure .  55 

5.  Convergence  of  column  buckling  loads .  58 

6.  Convergence  study  for  an  isotropic  plate  model .  59 

7.  Effect  of  element  aspect  ratio . 60 

8.  Deflections  of  an  equivalent  homogeneous  orthotropic  plate .  60 

9.  Deflections  of  a  cross-ply  [0/90]  laminate .  61 

10.  Convergence  of  unconstrained  and  constrained  critical  buckling  loads . 65 


ix 


INTENTIONALLY  LEFT  BLANK. 


1  Introduction 


Modern  structural  design  is  making  increasing  use  of  advanced  composite  material  systems 
to  improve  material  efficiency  by  minimizing  structural  weight  while  meeting  static  and 
dynamic  design  strength  criteria.  While  maximizing  structural  efficiency,  these  materials  ex¬ 
hibit  complex  elastic  and  failure  response  under  loading  such  as  anisotropy  and  multimode 
failure  mechanisms.  Of  particular  interest  is  the  susceptibility  of  these  material  systems  to 
impact  damage.  In  aerospace,  submersible,  and  armored  military  vehicle  design,  the  accurate 
prediction  of  residual  strength  due  to  impact  damage  is  important  for  assessing  survivability 
and  damage  tolerance  to  projected  damage  causing  impact  threats  in  the  intended  service 
environment.  With  an  increasing  need  to  reduce  production  costs  associated  with  experi¬ 
mental  testing,  numerical  simulation  is  being  emphasized  to  assess  evolving  design  concepts. 
A  robust  analytical  approach  is  thus  required  to  ensure  a  viable  structure  by  accurately 
assessing  the  effect  of  assumed  internal  damage  on  residual  strength. 

The  study  of  impact  involves  the  complex  dynamics  between  colliding  bodies  which  poten¬ 
tially  entails  high  energy  transfer  rates,  stress  wave  phenomena,  large  deformations,  material 
phase/property  change,  large  strains,  fragmentation,  fracture,  and  a  multitude  of  additional 
failure  mechanisms  many  of  which  are  material  specific.  Common  analytical  approaches  de¬ 
fine  high-  and  low-velocity  regimes  to  characterize  the  impact  event  and  in  which  to  apply 
various  simplifying  assumptions.  In  low  energy  impact,  damage  can  be  distributed  signifi¬ 
cantly  away  from  the  contact  zone  due  to  longer  contact  duration  in  which  large  fractions 
of  the  impact  energy  can  be  transferred  into  elastic  deflections  within  the  target.  In  high 
energy  regimes,  through-penetration,  explosive  spalling,  and  diffuse  fragmentation  can  char¬ 
acterize  the  material  damage.  If  the  ensuing  material  damage  has  compromised  structural 
integrity,  under  subsequent  loading,  damage  modes  will  propagate  until  total  structural  fail¬ 
ure  is  precipitated.  In  laminated,  fiber- reinforced  composites,  the  heterogeneous,  layered 
characteristics  of  the  material  system  exhibit  a  highly  complex  combination  of  microscopic 
and  macroscopic  damage  modes.  These  damage  modes  are  pronounced  in  low  velocity  impact 
events  and  interact  in  a  complex  manner  under  subsequent  loading.  Damage  modes  include 
material  failure  through  matrix  cracking,  fiber  breakage,  debonding  along  fiber-matrix  in¬ 
terfaces  and  ply  groups,  and  structural  failure  modes  associated  with  sublaminate  buckling 
due  to  delaminations.  Depending  on  the  geometric  complexity  of  the  problem  being  stud¬ 
ied,  various  approaches  have  been  used  to  model  the  problem.  Mathematical  models  have 
utilized  exact  formulations  based  on  elasticity  theory,  engineering  approximations  such  as 
beam,  plate  and  shell  assemblages,  and  numerical  approaches  such  as  finite  element  analysis 
[1-16].  For  a  complete  analysis,  both  the  local  and  global  response  of  the  component  to  ap¬ 
plied  loading  is  required.  This  response  includes  the  basic  displacement  and  internal  stress 
distribution  from  which  local  material  failure  modes  and  sublaminate  instabilities  may  be 
predicted  and  the  residual  strength  determined. 

Of  the  myriad  of  analytical  and  numerical  approaches  that  have  hitherto  been  used  to  pre¬ 
dict  residual  strength  in  composites,  these  analyses  have  generally  been  restricted  to  selected 
damage  modes  and  idealized  geometric  configurations.  Although  these  efforts  have  been  im¬ 
portant  in  elucidating  an  understanding  of  failure  phenomenon  in  composites,  they  have  not 
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permitted  a  direct  application  to  the  analysis  of  actual  structural  designs  due  to  inherent 
simplifications.  The  current  effort  has  been  directed  towards  incorporating  advanced  analyt¬ 
ical  representations  of  composite  damage/failure  phenomena  in  a  general  analysis  tool  which 
can  be  applied  to  actual  structural  geometries  and  realistic  damage  states. 

The  finite  element  method  provides  a  completely  generic  numerical  approach  to  potentially 
incorporate  all  salient  material  behavior  and  failure  modes  while  allowing  arbitrary  geome¬ 
tries,  applied  loads,  and  support  conditions  to  be  modeled  in  the  analysis  [17].  Prior  efforts 
directed  towards  predicting  residual  strength  using  a  finite  element-based  approach  have 
fallen  short  of  the  development  of  a  comprehensive  analysis  tool  [18-25].  Limitations  include 
restriction  to  two-dimensional  problems,  limited  selection  of  failure  criteria,  and  inability  to 
assess  both  material  and  structural  failure  modes  in  a  combined  progressive  analysis. 

This  report  details  a  numerical  finite  element-based  approach  for  predicting  residual  strength 
in  general  three-dimensional  structures  containing  internal  damage.  The  contribution  of  this 
effort  to  the  prediction  of  residual  strength  is  the  development  of  a  viable  solution  algorithm 
accounting  for  the  coupled  effects  of  nonlinear  material  response  and  both  material  and  insta¬ 
bility  failure  mechanisms  which  has  not  been  presented  in  the  literature.  While  emphasizing 
the  analysis  of  progressive  failure  in  laminated  composites,  the  support  of  a  layered  media 
permits  additional  material  systems  to  be  modeled,  such  as  sandwich-type  constructions, 
piecewise  linear  approximations  to  homogeneous  materials  with  varying  properties  along  a 
thickness  direction,  and  isotropic  materials  modeled  as  a  single  layer.  This  methodology  has 
been  incorporated  into  a  computer  program  designated  RESTRAN  (REsidual  STRength 
ANalysis)  for  determining  the  ultimate  strength  of  composite  structures  with  arbitrary  ge¬ 
ometric  configurations  and  general  loading  and  support  conditions.  The  flow  of  execution 
is  based  on  determining  a  sequence  of  scale  factors  to  the  initial  applied  loads  to  cause  a 
series  of  incremental  failures  leading  to  ultimate  collapse  of  the  structure.  The  sequence  of 
scale  factors  is  not,  in  general,  a  monotonically  increasing  function;  thus,  the  highest  load 
multiplier  attained  during  the  course  of  the  analysis  is  used  to  provide  a  measure  of  the  resid¬ 
ual  strength.  User-defined  subroutine  interfaces  have  been  incorporated  to  allow  alternative 
theoretical  approaches  to  be  applied  in  the  assessment  of  residual  strength,  which  thereby 
generalizes  RESTRAN  as  a  versatile  research  tool. 

The  following  sections  detail  the  finite  element  basis,  failure  modes,  damage  laws,  analy¬ 
sis  algorithms,  and  capabilities  of  the  RESTRAN  computer  program. 

2  Finite  Element  Formulation 

RESTRAN  has  been  developed  as  a  specialized  code  to  model  sequential  multimode  failure 
in  three-dimensional  composites.  The  development  of  an  extensive  library  of  solid  elements 
was  excluded  in  favor  of  a  single  eight-node  hexahedral  element  to  model  most  expected 
geometries  of  interest.  Higher-order  hexahedrals  were  also  avoided  in  the  present  effort  be¬ 
cause  of  the  expected  use  to  discretize  the  through-thickness  resolution  down  to  small  ply 
groups  or  even  local  regions  at  the  ply  level  such  that  a  local  linear  displacement  field  was 
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deemed  sufficient.  The  element  formulation  utilizes  the  hybrid  stress  technique  in  which 
both  stresses  and  displacements  are  assumed  as  independent  quantities  to  calculate  elastic 
stiffness  coefficients.  Incompatible  displacement  modes  are  introduced  to  complete  quadratic 
terms  associated  with  bending  and  condensed  into  the  element  formulation  to  enhance  el¬ 
ement  performance.  This  element  is  commonly  referred  to  as  the  Pian-Tong  Hexahedral 
[26].  To  model  geometric  nonlinear  response,  differential  stiffnesses  are  computed  based 
on  higher-order  strain  components  and  purely  displacement-based  shape  functions.  Within 
a  finite  element  framework,  structural  failure  through  sublaminate  buckling  is  accounted 
through  a  specialized  eigenanalysis  of  instability  modes  and  algorithmic  processing  of  mul¬ 
tiple  delaminations  in  which  contact  constraints  are  enforced.  The  exceptional  accuracy  of 
stress  recovery  provided  by  the  hybrid  formulation  yields  highly  accurate  differential  stiffness 
coefficients  which  lead  to  rapid  convergence  of  buckling  loads.  Material  degradation  due  to 
progressive  ply  failures  is  determined  using  selected  failure  criteria  and  specialized  ply  dam¬ 
age  laws  which  are  simulated  through  degradations  made  to  element  stiffness  coefficients. 

Basic  geometric  conventions,  elastic  and  differential  stiffness  formulations,  and  nonlinear 
material  and  buckling  solution  algorithms  are  presented  in  the  following  subsections. 

2.1  Coordinate  Systems  and  Nomenclature 

The  finite  element  formulation  in  RESTRAN  incorporates  several  conventions  pertaining  to 
coordinate  system,  laminate  description,  and  tensor  notation. 

2.1.1  Global  and  Local  Coordinates 

The  global  Cartesian  (x,y,z)  coordinate  system  is  used  for  describing  the  geometry  and  finite 
element  discretization  of  the  structural  model.  At  a  particular  element,  without  the  specifi¬ 
cation  of  a  user-defined  local  coordinate  system,  a  local  (x’,y’z’)  element  coordinate  system 
is  assumed,  which  is  offset  and  congruent  to  the  global  system  in  which  to  define  laminate 
orientation.  This  coordinate  system  is  depicted  in  Figure  1. 


Figure  1.  Global  and  local  coordinate  systems. 
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The  description  of  the  laminate  layup  assumes  that  ply  layers  are  parallel  to  the  local  x’y’- 
plane  and  that  the  stacking  sequence  is  a  function  of  the  local  z’-coordinate,  as  shown  in 
Figure  2(a).  The  principal  ply  directions  (1,2,3)  are  defined  by  the  offset  angle  from  the  x’- 
axis  and  are  referred  to  as  the  longitudinal,  transverse,  and  normal  coordinates,  respectively. 
These  quantities  are  depicted  in  Figures  2(b)  and  2(c). 


Figure  2(b).  Principal  ply  directions. 
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Figure  2(c).  Ply  orientation  angle  convention. 


Components  of  the  stress  and  strain  tensors  are  labeled  with  respect  to  the  global  (x,y,z) 
or  local  (x’,y’,z’)  coordinates,  as  shown  in  Figure  3.  When  transformed  to  principal  axes, 
the  principal  (1,2,3)  coordinate  designations  are  designated  as  longitudinal,  transverse,  and 
normal  directions,  respectively.  Stress  and  strain  vectors  are  ordered  as 

{&}  —  [Oxx,  @yyi  & zzi  &yzi  ®zn  &xy\  ^ -q 

{e}  —  [^ii,  ,  Cyy,  eZ2,  Cy2,  ezx,  Ciy] 

which  are  related  through  the  material  constitutive  law  expressed  in  general  as 


The  [C]  matrix  is  assumed  to  be  constituted  by  orthotropic  lamina  which  possess  the  x-y 
plane  as  a  plane  of  symmetry  and  are  described  by  nine  independent  material  moduli.  With 
arbitrary  orientation  of  the  constituent  layers,  the  assembled  material  stiffness  matrix  may 
assume  a  general  monoclinic  form  given  by 

C\2  C\z  0  0  Ci6 

C21  C22  C23  0  0  C26 

C31  C32  C33  0  0  Cze 

0  0  0  C44  C45  0 

0  0  0  C54  C55  0 

Cei  Cq2  Cqz  0  0  Cqq 
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2.1.2  Coordinate  Transformations 

A  local  material  coordinate  system  ( x ,  y,  z)  may  be  specified  for  each  element  in  RESTRAN 
to  map  element  stresses,  strains,  and  material  properties,  to  account  for  any  laminate  ori¬ 
entation  with  respect  to  the  local  element  (x’,y’,z’)  system.  These  Cartesian  coordinate 
systems  are  depicted  in  Figure  4. 


Figure  4.  Local  (x’,y’,z’)  and  ( x,y,z )  coordinate  systems. 


Transformation  matrices  are  defined  as 
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(4) 


[T£]  = 
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and  where  l,  m,  and  n  are  the  direction  cosines  between  the  two  Cartesian  systems  defined 
as 

(  x  )  T  k  mi  ni  1  f  x'  ) 

(7) 

J  UJ 

The  transformation  matrices  have  the  property  that 
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>  = 
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n2 

<  y’  > 

uj 
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IT,]-1  =  [T„f  [T,]-1 

which  yield  the  relations  for  stress  transformations 

=  [TJT 

(8) 

{*}  =  [Ta]{a’} 

K}  = 

[Te]T{<r} 

(9) 

and  relations  for  strain  transformations 

{*}  =  [T€]{6'} 

{e'}  = 

[T  „f{e} 

(10) 

and  finally,  for  material  transformations  as 

[C]  =  MC'][T/] 

[C]  = 

[T£t][C][T€] 

(11) 

2.2  Elastic  Stiffness  Formulation 

The  hybrid  stress  technique  is  utilized  to  form  elastic  stiffness  coefficients.  Details  on  the 
hybrid  method  may  be  found  in  References  [26-30].  The  structure  of  element  matrices  are 
defined  by  the  Hellinger-Reissner  functional  given  by 

IIr  =  J[{— 1/2)<ttS<7  +  crT(Lu9)  -  (Lt er)T \x\)dv  (12) 

where  a  is  the  assumed  stress  field,  S  is  the  material  compliance  matrix,  u9  and  ua  are  the 
assumed  compatible  and  incompatible  displacement  fields,  and  L  is  the  differential  operator 
relating  strains  to  displacements. 
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The  assumed  stresses  are  assumed  as 


a  =  P/3  (13) 

where  P  is  a  matrix  of  polynomial  terms,  and  f3  is  a  vector  of  undetermined  expansion 
coefficients.  The  displacement  field  is  assumed  over  the  element  domain  as 

u  =  u9  +  uA  (14) 

in  which  compatible  and  incompatible  displacement  components  are  given  by 

ug  =  Nq  (15) 

uA  =  MA  (16) 

where  N  and  M  are  compatible  and  incompatible  displacement  shape  functions,  respectively, 
q  are  nodal  displacements,  and  A  are  Lagrange  multipliers  which  enforce  internal  constraints. 
In  the  form  of  equation  (12),  performing  the  variation  with  respect  to  u\,  the  incompatible 
displacements  may  be  used  to  variationally  enforce  a  priori  the  field  equilibrium  conditions 
through  the  last  term  in  equation  (12)  which  requires  that 

6  J  (LT cr)T  u\dv  =  0  (17) 

or 

LTcr  =  0  (18) 

Applying  the  constraints  to  the  stress  modes  results  in  the  reduced  functional 

II#  =  J[(-1/2)<TTS(T  +  erT(Luq)]dv  (19) 

Substituting  equations  (13),  (15),  and  (16)  into  (12)  yields 

Hr  =  /[(— l/2)/3rPTSP/3  +  /3rPr(LN)q]du  (20) 

or 

UR  =  (— l/2)/3rH/3  +  /3rGq  (21) 

where 

H  =  j;PrSPdu  (22) 

G  =  fv  PT(LN)dv  =  fv  PTBdv  (23) 

Seeking  a  stationary  value  of  the  functional  by  taking  the  first  variation  with  respect  to  /3 
yields 

/3  =  H-1Gq  (24) 

Substituting  the  resulting  expression  for  (3  into  equation  (21),  the  variation  with  respect  to 
q  yields  the  element  stiffness  matrix  as 

K  =  GrH_1G  (25) 
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The  element  geometry  and  node  numbering  are  depicted  in  Figure  5. 


Figure  5.  Hexahedral  element  configuration. 


The  displacement  functions  ng  are  given  by 


uq  ] 

8  1 

U9  =  < 

,  =  Eo(1+  60(!  +  W)(l  +  CeC) < 

<= i  8 

Vi 

l  WQ  J 

{  wi  > 

(26) 


The  isoparametric  mapping  between  physical  and  natural  coordinates  is  given  by 


X  —  Oo  +  Gl£  +  02 7?  +  03^  +  &itV  +  +  C-eVC  "h  a7^VC 

y  =  bo  +  bit  +  b2T)  +  b3(  +  b4£rj  +  +  b6r)C  +  br&C  (27) 

2  =  c0  +  cit  +  c2r]  +  c3(  +  dtr]  +  c5tt  +  Cf>r)C  +  CrtvC 
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The  stress  field  is  initially  assumed  as  complete  quadratic  expansions  in  natural  or  tensor 
coordinates  given  by 
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where 


[Pol  =  [!,{,'),  C.ft.iK.Cf.fW.C2] 


(30) 


As  detailed  in  by  Pian  and  Tong  [26],  incompatible  displacement  functions  are  assumed  such 
that  the  polynomial  order  the  total  displacement  field  is  cubic  thereby  yielding  a  resultant 
quadratic  strain  field  of  the  same  order  as  the  assumed  stresses.  The  assumed  stress  modes 
are  then  subjected  to  the  constraint  given  by  equation  (18),  yielding  an  element  stress  field 
with  18  independent  stress  modes  given  by 
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(31) 


To  preserve  the  constant  stress  terms,  the  natural  stresses  are  mapped  to  physical  space 
through  a  contravariant  transformation  using  Jacobians  computed  at  the  element  centroid 


(32) 


The  stress  field  expressed  in  physical  or  Cartesian  coordinates  is  given  by 


where 


[P2}  = 


p = [  [Pi 

m  [p»i] 

(33) 
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The  eight-node  hexahedral  element  formulation  based  on  hybrid  element  theory  is  a  robust 
element  with  excellent  behavior  in  normal  stretching  and  bending  deformation  modes. 


2.3  Modeling  Material  Nonlinearity 

The  RESTRAN  analysis  program  supports  nonlinear  elastic  material  properties  and  nonlin¬ 
ear  inelastic  behavior  due  to  ply-level  failures.  The  combined  elastic  and  inelastic  nonlin¬ 
earities  are  accommodated  in  the  solution  algorithm  by  assuming  that  the  nonlinear  elastic 
response  is  independent  of  material  failure  and  maintains  the  same  differential  stress/strain 
relationship  at  each  point  during  unloading.  As  depicted  in  Figure  6,  any  reduction  in  mod¬ 
uli  due  to  the  occurrence  of  a  partial  failure  mode  is  assumed  to  result  in  a  stress/strain 
relation  self-similar  to  the  initial  nonlinear  elastic  curve.  These  two  forms  of  material  nonlin¬ 
earity  are  processed  sequentially.  The  primary  unknown  in  the  overall  RESTRAN  solution 
algorithm  is  the  scale  factor  applied  to  the  initial  set  of  external  loads  to  cause  the  next 
failure  event.  For  material  failure,  this  factor  is  determined  by  a  linear  extrapolation  using 
the  selected  failure  criteria  and  is  thus  a  function  of  the  current  stress-strain  state.  If  buck¬ 
ling  failure  is  included,  scale  factors  are  also  computed  to  cause  the  next  instability  failure. 
Thus,  convergence  to  the  desired  equilibrium  state  requires  stationarity  of  the  scaled  applied 
load  to  cause  next  failure  and  force  equilibrium  between  external  and  internal  loads  as  a 
function  of  the  nonlinear  stress-strain  relations.  Therefore,  the  algorithmic  implementation 
in  RESTRAN  first  converges  the  nonlinear  elastic  state  to  obtain  the  vector  of  external  loads 
to  cause  the  next  failure.  This  set  of  loads  is  then  used  to  compute  stresses  for  evaluation 
in  selected  failure  criteria.  A  typical  cascade  of  element  material  failures  is  subsequently 
accounted,  and  convergence  is  assumed  when  no  additional  ply  failures  occur  due  to  internal 
stress  redistribution  at  the  current  load  level. 

Assuming  small  strains,  at  some  applied  load  level,  F*,  the  stress/strain  relationship  may 
be  expressed  as 

<rk  =  Ck(e)ek  (37) 

where  Cfc(c)  is  the  strain  dependent  material  stiffness  matrix,  and  crk  and  ek  are  the  stress 
and  strain  vectors  at  the  kth  load  level.  For  greatest  generality,  the  particular  nonlinear 
stress-strain  relation  is  input  into  RESTRAN  via  a  user-defined  subroutine  as  described  in 
the  user’s  manual  [31]. 
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Figure  6.  Partial  failure  in  nonlinear  elastic  material. 


To  solve  the  nonlinear  system  of  equations,  an  iterative  scheme  is  employed  to  obtain  the 
equilibrium  state.  Because  convergence  is  sought  to  the  load  at  which  the  next  failure  event 
will  occur,  the  overall  nonlinearity  is  due  to  contributions  from  the  nonlinear  stress/strain 
relations,  the  extrapolations  of  the  selected  material  failure  criteria  which  may  be  quadratic, 
and  the  calculation  of  eigenvalues  associated  with  bifurcation  buckling  states.  The  potential 
volatility  of  the  solution  suggests  a  simple  iterative  approach  using  updated  secant  stiffnesses 
as  opposed  to  more  involved  tangent  matrix  methods  such  as  Newton-Rahpson  iteration  or 
matrix  update  schemes  such  as  the  BFGS  technique.  For  the  ktk  analysis  cycle,  iterations 
are  performed  such  that  the  sequence  of  loads,  F,,  converge  in  the  limit:  Fl5  F2,  F3,  — >,  Pk. 
In  the  current  implementation,  at  the  ith  iteration,  the  calculated  load  multiplier,  A*,  is  used 
to  scale  loads  converged  at  the  last  kth  analysis  cycle  to  obtain  a  new  set  of  external  loads 
as 

Fj  =  XiPk  (38) 

The  ith  approximation  to  the  total  displacements,  D,,  and  stress  dependent  secant  stiffness 
are  used  to  compute  the  internal  load,  Qi,  as 

Qi  =  KtDi  (39) 

The  residual  load  vector  measuring  the  difference  between  external  and  internal  loads,  Ri, 
is  obtained  from 

AR,  =  Fi  -  Qi  (40) 

The  increment  in  displacements  is  calculated  using  the  residual  as 

ADi+1  =  K-'ARi  (41) 
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and  the  new  estimate  for  total  displacements  satisfying  equilibrium  is  given  by 


Di+i  —  Dj  +  ADj+i  (42) 

In  the  course  of  determining  a  sequence  of  scale  factors  to  precipitate  the  next  failure,  the 
secant  moduli  act  to  stabilize  convergence  for  the  case  of  changing  external  loads  which  may 
increase  or  decrease  during  iterations  due  to  extapolations  made  from  nonlinear  failure  cri¬ 
teria.  This  scheme  is  depicted  in  Figure  7. 


This  nonlinear  solution  procedure  is  sufficient  for  weakly  to  moderately  nonlinear  systems. 
The  extrapolated  scale  factors  obtained  from  material  failure  criteria  and  the  eigenvalues  ob¬ 
tained  from  linear  stability  analysis  function  as  line  searches  towards  the  desired  minimum 
scale  factor  to  cause  the  next  failure  event. 

Iterations  are  continued  until  the  vector  norm  of  the  residual  load  vector  is  less  than  a 
prescribed  tolerance 

|AR*|  <Tol  (43) 

Because  nonlinearities  present  in  the  system  may  shift  the  target  applied  load  to  cause  the 
next  failure,  a  modification  of  the  current  scale  factor  is  made  to  damp  out  oscillations  and 
speed  convergence.  At  the  ith  iteration,  the  scale  factor,  A j,  is  formed  as  a  product  accounting 
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for  the  previous  estimate  for  the  next  failure  A;_i  and  the  scale  factor  obtained  directly  from 
the  selected  failure  criteria,  a;*.  This  relationship  is  expressed  as 


A  i  —  A  i—\QZi 


(44) 


The  estimate  for  cti  is  modified  using  a  biasing  factor,  pi,  to  speed  convergence.  This 
modification  is  given  by 


Aj  = 


Ai-i[l  +  (5  +  P){oii  -  1)]  for  cti  >  1 
Ai_i[ai  +  (|  -/?)(!-  a,-)]  for  a*  <  1 


(45) 


where  the  biasing  factor  is  formed  as  an  accumulated  sum  based  on  all  previous  estimates 
for  cti  as 


(46) 


and  subjected  to  the  restriction  that  0.1  <  Pi  <  0.9.  When  the  nonlinear  system  of  equations 
has  converged  to  an  equilibrim  state  corresponding  to  the  applied  load  Pk,  ply  failures  or 
sublaminate  instabilities  are  then  processed. 


2.4  Modeling  Geometric  Nonlinearity 


Geometric  instability  in  impact  damaged  laminated  structures  constitutes  an  important  and 
common  failure  mode  which  results  in  local  or  global  buckling  failure.  These  modes  become 
prevalent  when  the  damage  state  contains  significant  regions  of  reduced  stiffness  due  to  ma¬ 
terial  failure  or  extended  planes  of  delaminations  between  ply  groups.  Buckling  defines  an 
unstable  equilibrium  state  at  which  the  structure  cannot  elastically  resist  an  infinitesimal 
departure  from  the  current  configuration  under  the  applied  loads.  This  condition  is  charac¬ 
terized  by  the  change  in  the  total  potential  energy  associated  with  an  arbitrary  displacement, 
Sq,  from  equilibrium  is  an  extremum  or  equivalently,  that  the  total  potential  satisfies  the 
Trefftz  criterion  [32]  given  by 


n  =  u  +  v 


(47) 


and 


52n 

~W 


=  0 


(48) 


in  which  II  is  the  total  potential,  U  is  the  strain  energy,  V  is  the  potential  energy,  and  <f> 
is  the  generalized  displacement.  In  the  linear  theory  of  buckling,  a  solution  is  sought  at  the 
point  where  the  locus  of  equilibrium  states  bifurcate,  as  depicted  in  Figure  8. 
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Figure  8.  Bifurcation  of  equilibrium  states. 


In  order  to  approximate  the  nonlinear  geometric  buckling  behavior,  higher-order  strain  terms 
are  used  to  compute  the  differential  or  geometric  stiffness  matrix  Ka.  The  geometric  stiffness 
matrix  is  formulated  using  the  assumed  displacement  fields  in  (26)  and  is  given  by 


where 


and 


[Ka]  =  Jv[BNL)T[T][BNL]dv 


[r]  = 


[S]  0  o' 

0  [S]  0 

0  0  [S] 


O'xxO 

7~xy0 

Tzx  0 

1~xy0 

&yyO 

TyzO 

7~xz0 

TyzO 

&zzO 

(49) 


(50) 


(51) 
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The  B  hl  matrix  is  given  by 


where 


[Bnl, 


[G]  0  0 

0  [G]  0 

0  0  [G] 


"  N1)X 

0 

0 

n2jI  . 

••  nm- 

[G]  = 

Ni* 

0 

0 

n2„  . 

-  N  s. 

N1)Z 

0 

0 

..  n8jZ 

(52) 


(53) 


The  resulting  nonlinear  equilibrium  relationship  is  given  by 

(K  +  Kff)X  =  R  (54) 

Determining  the  load  level  at  which  bifurcation  occurs  requires  the  solution  to  the  following 
generalized  eigenvalue  problem 

(K  -  AK<r)X  =  0  (55) 

The  critical  buckling  load  is  then  given  by 

Kent  =  AR  (56) 

where  A  corresponds  to  the  lowest  positive  eigenvalue  of  the  global  system. 

For  large  systems  wherein  a  small  set  of  eigenvalues  are  sought,  algorithms  have  been  incor¬ 
porated  into  RESTRAN  based  on  subspace  iteration  [33-40].  The  basic  subspace  iteration 
utilizes  an  initial  set  of  iteration  vectors,  Xk,  which  are  selected  to  span  the  space  defined 
by  the  desired  eigenvectors,  #,  with  the  corresponding  eigenvalues,  A.  Iterations  are  then 
performed  to  converge  the  eigenvectors  and  eigenvalues  in  the  limit 

Ak  — *  A  and  Xk  — >•  d?  as  k  — oo  (57) 

The  subspace  mapping  is  defined  by 

KX*+1  =  KaXk  (58) 

The  mapping  or  iteration  vectors  are  selected  in  accordance  with  procedures  outlined  by 
Bathe  [33].  The  selection  is  specified  such  that  the  first  column  in  X  is  taken  as  the  diagonal 
in  Kq,  while  remaining  columns  are  assigned  a  unit  value  in  the  row  position  corresponding 
to  min(kii/kau).  The  global  elastic  and  differential  stiffness  matrices  are  then  mapped  into 
the  subspace  as 

Kfc+1  =  XTk+lKXk+1  (59) 

Kak+1  =  Xl+lK„Xk+1  (60) 

resulting  in  the  reduced  eigenvalue  problem  given  by 

Kfc+iQfc+i  =  K(7fc+1  Qfc+iAfc+i  (61) 
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where  A  is  a  diagonal  matrix  of  the  system  eigenvalues.  Because  the  initial  iteration  vectors 
are  usually  approximate,  an  improvement  to  the  system  eigenvectors  is  given  by 

X-k+i  =  Xfc+1Q^+a  (62) 


and  an  iteration  of  the  subspace  is  continued  until  the  desired  set  of  eigenvalues  at  the 
(k  +  l)st  iteration  reach  a  stationary  value  defined  by  a  tolerance  as 


|A*+1  -  Af| 
|Af+1| 


<  tol 


i  =  1,2,3  ,...,n 


(63) 


During  iterations,  a  check  is  made  to  ensure  that  the  smallest  eigenvalues  are  being  obtained 
by  using  the  Sturm  sequence  property  of  the  characteristic  polynomials  of  eigensystems.  If 
a  shift  /i  is  introduced  slightly  greater  than  the  current  estimate  of  Xn,  the  factorization 


Ks  =  K  -  fiKa-  =  LDLt 


(64) 


yields  a  diagonal  matrix  D,  in  which  the  number  of  negative  diagonals  is  equal  to  the  number 
of  eigenvalues  smaller  than  p. 


An  additional  eigenvalue  extraction  technique  has  been  implemented  in  RESTRAN,  which 
is  based  on  the  Arnoldi  method  and  is  a  powerful  extension  of  subspace  iteration  [34,35]. 
After  reducing  the  generalized  eigenvalue  problem  definition  to  a  standard  format  through 
a  Cholesky  reduction  of  K  =  LTL,  pre-  and  post-multiplying  such  that 


A  =  L~XKCL  (65) 

yields  a  symmetric  A  matrix.  Arnoldi’s  method  for  finding  a  few  eigenvalues  of  A  proceeds 
as  follows:  given  an  initial  vector  xx  with  unit  norm,  at  each  step  m  construct  an  orthonormal 
basis  Xm  =  [xi,x2,  ...,xm]  for  the  Kyrlov  subspace  Km  spanned  by  [xi,  Axx, Am-1xx] 
by  computing  w  =  Axml  and  orthonormalizing  w  with  respect  to  xx,  x2, ...,  xm_x  to  obtain 
xm.  The  matrix  H  =  XTAX  is  upper  Hessenberg,  and  its  eigenvalues  provide  approxima¬ 
tions  to  m  eigenvalues  of  A.  The  main  steps  in  the  basic  Arnoldi  algorithm  for  finding  A j 
may  be  summarized  as  follows: 

1.  Initialization-.  Choose  the  number  of  steps  m  and  an  initial  vector  xx  with  unit  norm. 

2.  Arnoldi  steps  for  j  =  1,  2,  ...  ,  m: 

(a)  wj  =  Axj 

(b)  hy  =  x^Wj,  i  =  1,  2,  ...,j. 

(c)  sj  =  wj  -  £j=x  Xihy 

(d)  hj+Xj  =  1 1  ® j  1 1 2  >  xj+i  =  Sj/hj+ij- 


Set  X  =  [xx,x2,...,xm[ 


3.  Eignevalue  computation :  Reduce  the  upper  Hessenberg  matrix  H  =  {hy}  to  real  Schur 
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form  T  =  ZTHZ,  where  T  is  a  block  triangular  matrix  with  the  eigenvalues,  A;,  ordered  in 
descending  order  of  their  absolute  values  along  the  diagonal.  Set  X  =  XZ. 

4.  Convergence  test:  If  the  first  column  x,  of  X  satisfies  the  convergence  criteria  for  the 
residual 


reSi 


||(AX-XT)i||a 

l|A||a 


then  accept  Ai ;  otherwise,  return  to  step  2  and  repeat. 


(66) 


In  a  complex  structure,  the  various  load  paths  will  yield  positive  and  negative  eigenvalues 
corresponding  to  load  multipliers  to  cause  a  buckling  failure.  For  the  automated  interpreta¬ 
tion  of  instability  failure,  if  all  eigenvalues  are  negative  or  complex,  no  reversal  of  input  loads 
is  considered;  it  is  assumed  that  either  the  initial  set  of  applied  loads  preclude  a  buckling 
response  or  that  element  failures  have  progressed  to  a  terminal  failure  state.  In  the  event  of 
severe  numerical  difficulty  in  obtaining  eigenvalues  due  to  widespread  material  degradation, 
the  stability  analysis  is  bypassed,  and  any  further  failure  events  prior  to  prediction  of  catas¬ 
trophic  collapse  are  accounted  through  material  failure  modes. 


3  Laminated  Material  Representation 

RESTRAN  is  formulated  to  represent  layered  materials  such  as  laminated  composites.  In 
such  a  representation,  the  material  is  assumed  to  consist  of  an  assemblage  of  orthotropic 
laminae  or  plies.  The  discretization  of  the  model  into  finite  elements  associates  the  element 
properties  with  a  single  ply  or  with  a  number  of  plies  forming  a  local  sublaminate  group. 


3.1  Effective  Material  Moduli 

An  effective  approximation  to  the  the  mechanical  properties  of  an  assembled  sublaminate  is 
obtained  by  averaging  ply  properties.  For  thin,  plate-like  structures  under  long-wavelength 
loading,  classical  lamination  theory  (CLT)  has  been  adequate  for  obtaining  overall  mechan¬ 
ical  response  [41,42],  For  thick  laminates,  however,  more  accurate  formulations  for  effective 
three-dimensional  properties  have  been  developed  [43-45]  of  which  the  formulation  of  Sun 
and  Liao  [45]  has  been  incorporated  into  RESTRAN. 


The  sublaminate  constitutive  relationship  is  assumed  as 
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where  the  [C]  matrix  is  assumed  to  be  constituted  by  lamina  which  may  be  generally  or¬ 
thotropic  and  described  by  nine  independent  material  moduli.  In  general,  the  assembled 
material  stiffness  matrix  may  possess  only  symmetry  about  z  =  0  and  assume  an  anisotropic 
monoclinic  form  given  by 
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The  effective  three-dimensional  components  0^  are  given  by 


(68) 


Cu  =  E  VkQn  +  E(0?s  -  Au)«4((51,  -  %)/<& 

fc=l  k= 2 


(69) 


Ca  =  E  +  E(Qiz  -  Ai3 )<*(<&  -  <&)/<3» 

k=  1  A;=2 
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where 


(80) 


(81) 

(82) 

(83) 

The  transformed  stiffnesses  for  the  kth  ply  are  given  by 

Qn  =  Q\lcosA9k  +  Ql2sinA9k  +  2Q\2sin29kcos29k  +  4Qg6cos29k  sin29k 
Qk22  =  QknsinA9k +  Qk22cos49k +  2(Qkn  +  2Qk66)sin29kcos29k 

Qls  =  Q\z 

Qu  =  Qki4cos29k  +  Qk5sin29k 
0s5  =  Qk55cos29k  +  Qk44sin29k 

Qle  =  Qkee  +  (Qku  +  Qk22-2Qk2-4Qk6)sin29kcos29k  (84) 

Q12  =  Qk12  +  (Qkn  +  Qk22-2Qk2-4Qke)sin29kcos29k 
Qk 3  =  Qkl3cos29k  +  Qk3sin29k 
Q\ 3  =  Q\3sin29k  +  Qk23cos29k 

Qk6  =  (<3*lCos20fc  -  Q22sin29k  -  cos29k(Qk2  +  2 Qk66))cos9ksin9k 
Qk 6  =  (Qknsin29k -Q$2cos29k +  cos29k(Qk2  +  2Qk66))sin9kcos9k 
Qle  =  (Q\3-Qk23)sin9kcos9k 
Qk5  =  ( Qk5-Qk4)sin9kcos9k 

in  which 

Q"u  =  E*(  1  -44  )/(ElEl&) 

Qi 2  =  B2‘(4  +  44)/(£fB3‘fi‘) 
of3  =  £‘(4  +  44)/(B‘B‘n‘) 

<&  =  £*(1  -  44i)/(£f£3‘fi*) 

=  £»(>4>  +  ^i4)/(£f-S2n*)  (85) 

<4  =  s‘(i  -  44)/(£f£2tn‘) 

Q/c  s~ik 

44  “*  ^23 

<&  =  c* 

C?66  =  ^12 


4  =  4(^fc /£?) 
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(86) 


4  =  ^(El/E*) 

4  =  4  (£*/£?) 

^  =  (1  ~  ^12^21  ~  ^23^32  —  ^13^31  —  ‘^l/21l/32l'l3,) 

Eij ,  (7j;,  i/i;  and  0  are  the  Young’s  normal  and  shear  moduli,  Poisson  ratios  and  fiber  orien¬ 
tation  angle,  respectively. 

In  the  case  of  a  composite  material  system  with  a  balanced  [±0]  layup,  the  effective  material 
properties  reduce  to  those  obtained  using  classical  three-dimensional  lamination  theory  in 
which  the  material  stiffnesses  coefficients,  C y,  are  obtained  simply  as 

c a  =  i  E  ht Q»  (87) 

■"  k= 1 

where  H  is  the  local  sublaminate  thickness,  h*  are  the  individual  ply  thicknesses  and  N  is 
the  total  number  of  plies  in  the  sublaminate. 

The  assumption  of  a  layered  medium  permits  the  modeling  of  a  large  class  of  material 
systems:  laminated  composites,  sandwich-type  constructions,  piecewise  linear  approxima¬ 
tions  to  homogeneous  materials  with  varying  properties,  and  isotropic  materials  represented 
as  a  single  ply. 


3.2  Ply-Level  Stress  and  Strain  Recovery 

The  effective  material  constitutive  relationship  developed  by  Sun  and  Liao  [45]  is  based  on 
the  assumption  of  long  wavelength  loading  and  that  the  thickness  of  the  sublaminate  is  small 
compared  to  the  total  laminate  thickness.  This  leads  to  the  consequence  that  in  the  recovery 
of  average  or  effective  stresses  and  strains  at  the  element  level  using 

{1}  =  |B]{u}  {<r}[C]  =  {«}  (88) 


where  [B]  is  the  strain-displacement  matrix  and  {u}  is  the  displacement  vector,  the  inplane 
strains  and  transverse  normal  stresses  are  assumed  constant.  This  condition  may  be  stated 
as 

^XX 

®zz 
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K}  =  (C‘]{e‘) 


(90) 


the  following  equation  for  the  out-of-plane  stresses  results  in 
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and  the  unknown  out-of-plane  strain  components  may  then  be  calculated  from 


€zz  1 

Czz 

C34 

C35  ' 

-~1 

r 

f  1 

'Cn 

£ 

to 

C36 

\  €zz  1 

> 

1 

eyz  i 

>  rz 

C43 

C44 

C45 

< 

1  „  1 

|  aVz  ! 

►  — 

c4l 

C42 

C46 

< 

>  > 

j 

k 

.  C53 

C54 

C55 

k 

1 

[  Oxz  J 

k 

.  C51 

C52 

C$6  _ 

k 

1  ««  J 

k > 

(92) 


With  the  complete  strain  vector  for  the  kth  ply  layer  determined,  the  remaining  stress  com¬ 
ponents  are  obtained  using  equation  (90). 


3.3  Nonlinear  Elastic  Moduli 

As  shown  in  Figure  9,  nonlinear  elastic  materials  are  those  which  exhibit  a  nonlinear  stress/ 
strain  response  under  applied  loads.  Many  nonlinear  materials  exhibit  an  elastic  response 
below  the  fracture  or  yield  limit  in  which  deflections  follow  the  same  load-displacment  curve 
without  permanent  deformation  under  unloading.  The  material  moduli  for  these  materials 
may  therefore  be  expressed  as  continuous  functions  of  the  stress  or  strain  state  as 

Cjj  =  Cjj(<T*,e*)  (93) 

where  a*  and  e*  represent  individual  stress  or  strain  components  or  combined  measures  that 
represent  unit  potential  energy  or  strain  energy  density. 


Figure  9.  Nonlinear  elastic  stress-strain  relation. 
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For  greatest  generality  in  defining  nonlinear  properties,  nonlinear  stress/strain  relations  are 
input  as  secant  moduli  from  a  required  user-defined  subroutine.  This  subroutine  may  be  writ¬ 
ten  in  FORTRAN  or  C,  and  then  compiled  to  link  with  the  main  RESTRAN  executable. 
Through  a  standardized  arqument  interface,  RESTRAN  provides  stresses  and  strains  inter¬ 
polated  to  each  ply  and  at  element  integration  or  Guass  points.  The  user-defined  routine 
must  calculate  reduction  factors,  which  are  used  to  obtain  current  measures  of  the  secant 
material  moduli  depicted  in  Figure  10.  The  secant  moduli  are  computed  as 


Ef  = 

El  =  t 
El  =  >kEl 
<%  = 

<?31  =  ’/’5(?31 

G?2  = 


(94) 


A  full  description  of  this  subroutine  is  contained  in  the  RESTRAN  User’s  Manual  [31]. 


Figure  10.  Secant  moduli  depiction. 
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4  Material  Failure  Modes 


Composite  materials  exhibit  a  wide  range  of  failure  modes  including  microscopic  fiber  break¬ 
age,  shear-out  failure,  kink  band  formation,  and  microbuckling,  matrix  tensile  cracking  and 
compressive  failure,  and  macroscopic  multilayered  delamination  growth.  The  RESTRAN 
analysis  code  models  various  material  failure  modes  at  the  element  ply  level  using  selected 
failure  criteria  and  associated  damage  laws.  The  current  implementation  of  material  failure 
prediction  is  restricted  to  criteria  defined  at  a  point.  Nonlocal  failure  criteria  based  on  av¬ 
eraged  stress  or  strain  values  necessarily  require  the  definition  of  characteristic  integration 
paths  which  would  extend  over  many  element  domains  and  is  difficult  to  automate  in  a 
progressive  failure  analysis  without  introducing  some  degree  of  problem  specificity.  In  ad¬ 
dition,  failure  modes  associated  with  delamination  extension  are  currently  not  accounted  for. 

In  order  to  maximize  versatility  by  allowing  specialized  tailoring  of  the  analysis,  RESTRAN 
incorporates  user-defined  subroutine  interfaces  to  allow  any  material  failure  criterion  and 
damage  law  to  be  applied.  The  format  for  these  subroutines  are  discussed  in  detail  in  the 
RESTRAN  User’s  Manual  [31]. 


4.1  Material  Failure  Criteria 

A  primary  consideration  in  the  analysis  of  residual  strength  in  composite  structures  is  the 
accurate  prediction  of  ply  failure  due  to  a  specific  state  of  stress.  In  composite  materials,  the 
micromechanical  complexity  of  the  ply  material  makes  the  analysis  of  local  failure  difficult. 
Under  applied  external  loading,  the  local  stresses  contribute  to  a  large  number  of  identifiable 
failure  interactions  and  damage  modes,  any  one  of  which  can  become  a  critical  £weak-link’  in 
the  prediction  of  ply  failure.  A  large  number  of  criteria  have  been  derived  to  relate  internal 
stresses  and  experimental  measures  of  material  strength  to  the  onset  of  failure  [46-58].  All 
criteria  may  be  classed  as  degenerate  cases  of  the  general  tensor  polynomial  failure  criterion 
which  has  been  developed  with  the  general  form  given  by 

FI  =  Flol  +  FijOiOj  +  FijkOi<jjOk  +  ...  (95) 

where  cr,  are  stress  tensor  components  in  principal  material  directions,  Fu  Fij,  Fijk  are 
components  of  the  strength  tensors,  and  FI  is  the  failure  index  which  predicts  failure  when 
FI  >  1.  For  practical  applications,  the  terms  in  equation  (95)  higher  than  quadratic  are  usu¬ 
ally  set  to  zero  due  to  the  diminishing  returns  of  including  higher-order  polynomial  terms 
in  predicting  failure  together  with  the  increased  cost  of  experimentally  determining  these 
higher-order  strength  components. 

Developed  failure  criteria  exhibit  a  wide  range  of  predictive  capability.  The  simplest  utilize 
individual  stress  or  strain  component  relations  in  which  each  component  is  assumed  to  act 
independently  while  other  criteria  account  for  coupling  and  interactions  in  the  stress/strain 
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state.  In  addition,  many  criteria  simply  predict  an  unidentified  failure  state,  while  others 
differentiate  between  competing  failure  modes.  For  use  in  failure  criteria,  principal  ten¬ 
sor  components  are  computed  for  each  layer  in  laminated  composites  along  longitudinal, 
transverse,  and  normal  ply  coordinates.  For  laminates  exhibiting  monoclinic,  orthotropic, 
or  transversely  isotropic  material  properties,  global  stresses  in  the  local  (x’,y’,z’)  coordinate 
system  are  first  mapped  into  the  local  material  ( x,y,z )  coordinate  system  if  the  laminate 
is  arbitrarily  oriented  in  space.  Next,  principal  ply  stresses  in  the  lamina  (1,2,3)  coordi¬ 
nate  system  are  obtained  using  direction  cosines  which  are  determined  soley  by  the  fiber 
orientation  angle.  This  transformation  is  given  by 
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Principal  strains  are  obtained  similarity  as 
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The  direction  cosines  are  given  by 


l\  =  cos  3  mi  —  sin/3 

k  =  —sin/3  m2  =  cos/3 


(96) 


(97) 


(98) 


where  /3  is  the  fiber  orientation  angle. 

For  homogeneous  isotropic  materials,  the  principal  stresses  are  determining  by  projecting 
the  six  components  of  stress  onto  a  plane  where  the  shear  stress  components  vanish.  This 
leads  to  the  relation 

(P  Txy  TXZ 

~Txy  ip  ~  &y)  ~'ryz 
Txz  Tyz  (<7  Oz ) 

from  which  principal  stresses  and  principal  planes  may  be  determined.  An  additional  geo¬ 
metric  constraint  relates  the  direction  cosines  l,m,n  as 

l2  +  m2  +  n2  =  l  (100) 
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For  a  nontrivial  solution,  the  determinant  of  the  system  given  in  equation  (99)  must  vanish, 
which  leads  to  a  cubic  equation  for  the  principal  stresses  given  by 

O3  “b  A2(72  +  AjO"  -f-  Aq  =  0  (101 ) 
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To  obtain  the  solution,  the  following  quantities  are  defined 


r  —  -(AiA2  —  3A0)  —  —  A23 

*i  =  [r  +  (g3  +  r2)^ 

«2  =  [r  ~  (<J3  +  f2)^}  2 


(102) 


(103) 


(104) 


in  which,  for  stress  and  strain  tensor  transformations,  the  quantity  q3  +  r2  >  0  thus  guaran¬ 
teeing  that  all  roots  are  real.  The  principal  normal  stresses  are  then  given  by 

°i  =  (*i  +  s2)  -  y 

02  —  —  -(S1  +  S2) — ^  + -^—(si  —  s2)  (105) 

1/  ,  \  0,2  i\/ 3,  s 

°z  —  “2(si  +  52)--2 - 2~(3i  -  «2) 

Principal  shear  stresses  may  be  expressed  in  terms  of  the  principal  normal  stresses  as 

(74  =  |<72  —  <73 1 

^5  =  kl  -  <73 1  (106) 

0&  =  |(7x  —  02 1 


A  variety  of  failure  criteria  have  been  incorporated  into  the  RESTRAN  analysis  code  to 
predict  ply  failures  under  external  applied  loads  and  are  described  in  the  following  subsec¬ 
tions.  These  criteria  have  been  used  to  form  the  central  material  failure  analysis  capability 
in  other  state-of-the-art  analysis  codes  [59].  Strains  are  recovered  at  the  Gauss  points  and 
interpolated  to  each  ply  at  the  centroid  from  which  layer  stresses  are  obtained.  It  is  at  these 
locations  that  the  available  failure  criteria  are  applied.  For  greatest  generality,  a  user-defined 
subroutine  may  be  used  to  incorporate  any  specific  failure  criterion.  Both  interpolated  and 
Gauss  point  stresses  and  strains  are  passed  into  the  routine.  This  subroutine  is  described  in 
the  RESTRAN  User’s  Manual  [31]. 
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4.1.1  Maximum  Stress  Criterion 


Of  the  various  anisotropic  failure  criteria,  the  maximum  stress  criterion  combines  early  inves¬ 
tigations  into  failure  mechanisms  due  to  Rankine  and  Tresca.  The  failure  strength,  af,  may 
be  based  on  fracture  strength,  yield  strength,  proportional  limit,  endurance  limit,  maximum 
working  stress  or  some  other  scalar  parameter,  depending  on  the  expected  failure  mode  or 
design  criterion.  The  failure  modes  predicted  by  this  criterion  are  each  dependent  on  only 
one  stress  component  and  are  summarized  in  Table  1.  In  general,  for  ductile  materials, 
the  maximum  normal  stress  criterion  yields  poor  failure  predictions,  while  the  maximum 
shearing  stress  predictions  produce  good  results.  For  brittle  materials,  the  maximum  stress 
criterion  generally  produces  good  predictions  of  failure  strength.  In  applying  this  criterion, 
ultimate  ply  failure  is  predicted  when  any  one  of  the  following  conditions  are  satisfied 
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where  <t2 ,  a2,  and  a3  are  the  principal  normal  stress  components  and  a4,  as,  and  a6  are  the 
principal  shear  stress  components.  The  strength  measures  are  given  by  Xt,  Xc,  Yt,  Yc,  Zt, 
and  Zc,  which  are  the  normal  tensile  and  compressive  strengths  in  the  principal  1,  2  and  3 
directions,  respectively,  and  R,  S,  and  T  are  the  shear  strengths  defined  in  the  23,  13,  and 
12  planes,  respectively.  Expressed  in  tensor  polynomial  form,  the  maximum  stress  criterion 
may  be  written  as 

(a1-Xt)(a1+Xc)(a2-Yt)(a2  +  Yc)(a3-Zt)(a3  +  Zc)(\a4\-R)(\as\-S)(\a6\~T)  =  0  (108) 


Table  1.  Maximum  stress  criterion  failure  modes. 


Condition 

Failure  Mode 

o\  >  Xt  and  o\  >  0 
0i  <  Xc  and  0i  <  0 
02  >Yt  and  <72  >  0 
02  <YC  and  <j2  <  0 
03  >  Zt  and  a3  >  0 
a3  <  Zc  and  a3  <  0 

kl  >  1 
kl  >  1 
kl  >  1 

Fiber  failure  in  tension 

Fiber  failure  in  compression 
Matrix  failure  in  tension 
Matrix  failure  in  compression 
Matrix  failure  in  tension 
Matrix  failure  in  compression 
Interlaminar  shear  failure 
Interlaminar  shear  failure 
Inplane  shear  failure 

4.1.2  Maximum  Strain  Criterion 

The  maximum  strain  criterion  is  of  identical  form  as  the  above  with  stresses  replaced  by  cor¬ 
responding  strains  and  strength  measures  replaced  by  associated  ultimate  strain  measures. 
This  criterion  is  founded  on  early  developments  made  by  St.  Venant.  As  with  the  maximum 
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stress  criterion,  a  number  of  failure  modes  may  be  identified  as  listed  in  Table  2.  Ultimate 
ply  failure  is  predicted  when  any  one  of  the  following  conditions  are  satisfied 
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> 
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where  ex,  e 2,  and  e3  are  the  principal  normal  strain  components  and  e 4 ,  €5,  and  are  the 
principal  shear  strain  components.  The  ultimate  strain  measures  are  given  by  Xt,  Xc,  Yt,  Yc, 
Zt,  and  Zc  which  are  the  normal  tensile  and  compressive  strains  in  the  principal  1,2,  and  3 
directions,  respectively,  and  R,  S,  and  T  are  the  shear  strains  in  the  23,  13,  and  12  planes, 
respectively. 

Table  2.  Maximum  strain  criterion  failure  modes. 


Condition 

Failure  Mode 

ei  >  Xt  and  ex  >  0 
ex  <  Xc  and  ex  <  0 
e2  >  Yt  and  e2  >  0 
e2  <  Yc  and  e2  <  0 
e3  >  Zt  and  e3  >  0 
€3  <  Zc  and  e3  <  0 
|e4|  >  1 
ie5  i  >  1 
|ee|  >  1 

Fiber  failure  in  tension 

Fiber  failure  in  compression 
Matrix  failure  in  tension 
Matrix  failure  in  compression 
Matrix  failure  in  tension 
Matrix  failure  in  compression 
Interlaminar  shear  failure 
Interlaminar  shear  failure 
Inplane  shear  failure 

4.1.3  Beltrami  Criterion 

The  Beltrami  failure  criterion  assumes  material  isotropy  and  is  based  on  comparing  the  total 
strain  energy  per  unit  volume  of  a  multiaxial  stress  state  with  the  strain  energy  produced 
by  a  uniaxial  test  at  failure.  This  criterion  is  given  by 

1 

oc  =  [<7i2  +  C22  4-  cr32  -  2i/(<7i cr2  +  a203  +  CT3C1)] 2  (110) 

and  failure  is  predicted  when  the  combined  stress  measure  equals  or  exceeds  the  critical 
stress,  cry,  to  cause  failure  in  a  uniaxial  test. 

oc  >  oy  (111) 


4.1.4  Von  Mises  Criterion 

The  Von  Mises  failure  criterion  assumes  material  isotropy  and  is  based  on  a  decomposition 
of  the  strain  energy  density  into  volumetric  or  dilatational  energy  and  distortional  energy. 
This  combined  stress  measure  assumes  that  failure  in  a  multiaxial  stress  state  is  due  solely 
to  the  distortional  energy  of  the  system.  The  resulting  failure  criterion  may  be  expressed 
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using  principal  stresses  as 


°c  —  ~  0h)2  +  (°2  —  C3)2  +  (a3  —  0i)2]2  (112) 

wherein  failure  is  predicted  when  the  combined  stress  measure  exceeds  the  ultimate  yield 
strength  of  the  material  or 

crc  >  ay  (113) 


4.1.5  Hoffman  Criterion 

The  three-dimensional  Hoffman  criterion  [52,57]  is  given  by 

FI  =  F\(J\  +  F2&2  +  T3O3  +  £n<7j  +  Fi2<7i<72  +  F\3<7i<J3  +  i'22cr2 

F +  i?33<73  +  -^44^4  +  £5505  +  F66al 
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For  two-dimensional  plane  problems,  assuming  <73  =  <713  =  023  =  0  and  due  to  symmetry 
about  the  2,3  axis,  Y  =  Z,  the  Hoffman  criterion  reduces  to 


FI=( 


Xt  X , 


1  ,  .1  1  , 
m  +  (y  -  y>G 2  + 


01 


XX 


0^2 2  £12  _ 

Ttyc  +  T2  xtxc 


(116) 


4.1.6  Hill  Criterion 

The  three-dimensional  Hill  criterion  [50]  is  based  on  exclusively  quadratic  terms  given  by 
FI  =  Fiio\  +  F 12 O'] (72  +  F 130^10^3  +  ■F’220’2  +  F23 <72(73  +  -£33(73  +  £44(74  +  £55(75  +  F§§g\  (117) 


where 
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in  which  the  value  for  the  normal  strengths  X,  Y,  and  Z  are  taken  as  compressive  or  tensile, 
depending  on  the  sign  of  the  corresponding  normal  stresses.  For  two-dimensional  plane 
problems,  assuming  o3  =  ai3  =  023  =  0  and  due  to  symmetry  about  the  2,3  axis,  Y  =  Z, 
the  Hoffman  criterion  reduces  to 


FI  —  Fn(j\  +  F12O1G2  +  F22O2  +  F3  e<7g 


(119) 


4.1.7  Tsai-Wu  Criterion 

The  three-dimensional  Tsai-Wu  criterion  [53]  is  given  by 

FI  =  F\(7i  +  F2U2  +  F3<73  +  Fu<jf  +  F\2<J\02  +  F\3(J\(T3  +  F22G2  T 
F2$02&3  +  F33g\  +  F 44  +  F33o\  +  F33(Tq 


(120) 
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Fiz  =  --/(y/XtXeZtZc) 

F23  =  ~/(y/YtYeZtZe) 

For  two-dimensional  plane  problems,  assuming  a  slightly  different  form  for  the  interaction 
term  F12,  the  Tsai-Wu  criterion  reduces  to 


7-.r  ,  1  1  .  ,1  1  .  <7i2  022  <J\2 
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4.1.8  Christensen  Criterion 

The  Christensen  failure  criterion  differentiates  failure  into  three  possible  inodes  [48].  These 
are  fiber  tension,  fiber  compression,  and  matrix  failure.  Fiber  failure  in  tension  is  predicted 
when 


(123) 


(124) 


(125) 


A,  =  a(l  —  21/12)  /  {P2  Ei) 

A2  =  o(l  —  V21  ~  v2z)  j  Ef) 

A3  =  2(1  +  i/12)2/(3^2£'i2) 

An  =  2[(1  +  v2\  4-  ^2i2)  +  ^23(1  —  ^21)  +  l/2s2]/(3/32E22)  (126) 

Ab  =  2[(2z/2i  +  ^23) (1  +  z/i2)]/(3^2£'i£,2) 

As  =  2[(— 1  +  2^21  +  2i/2i2)  —  21/23(2  +  U21)  —  v2z2]/ (3/32E22) 

Ay  =  1/ (2/?2(?232) 

As  =  1/(2/92Gi22) 

and  where  the  constants  a  and  (3  are  determined  experimentally  for  the  particular  material 
system  being  used. 


Xt 


1 01  —  ^1202  —  ^1303 
Fiber  failure  in  compression  is  predicted  when 


<  1  for  (oi  -  1/1202  -  ^1303)  >  0 


<1  for  (ai  -  uua2  -  1/1303)  <  0 

1 01  —  ^1202  —  ^1303  I 

Matrix  failure  is  predicted  when  the  following  condition  is  satisfied 

A\0'  1  +  -/4.2 (cr 2  +  03)  +  ^30i2  +  A^ipf2  +  032)  T 
44501  (02  +  03)  +  ^60203  +  AyO^  +  -4s(052  +  Cq2)  >  1 


where  the  constants  ^4i  — >  As  are  given  by 


4.1.9  Feng  Criterion 

The  Feng  failure  criterion  is  based  on  strain  invariants  and  predicts  fiber  breakage  and  matrix 
cracking  failure  modes  [47] .  Matrix  failure  is  predicted  when 


A\J\  +  A2J\‘  +  A3J2  A  1 
and  fiber  failure  is  predicted  when 


(127) 
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A4J5  +  A§J^  +  ^46=^4  >  1  (128) 

In  the  above  equations,  4li  — >  A$  are  experimentally  determined  constants,  and  J;  are  strain 
invariants  given  by 

J\  —  Ci  +  £2  +  63 

J2  =  g[(fl  —  C2)2  +  (^3  —  C2)2  +  (€1  —  C3)2  +  C42  +  C52  +  ^62] 

J4  =  e42  4-  e52  (129) 

Jb  =  €x 


4.1.10  Modified  Hashin  Criterion 

The  modified  Hashin  failure  criterion  [49]  differentiates  between  fiber  and  matrix  failure 
modes.  Fiber  failure  is  predicted  on  the  basis  of  whether  the  fiber  or  longitudinal  stresses 
exceed  the  maximum  fiber  directional  strength.  Thus,  fiber  failure  in  tension  is  predicted 
when 


tt  >  1  for  ox  >  0  (130) 

■X-t 

and  fiber  failure  in  compression  is  predicted  when 

-^7-  >  1  for  0i  <  0  (131) 

Matrix  failure  is  predicted  using  maximum  normal  and  transverse  stress  components  in  the 
2-3  plane  defined  by 

Onn  =  (--^— )  +  (~  2  —)cos(2P)  +  a23sin(2/3) 

ont  =  2  —)sin{2j5)  +  o~23cos(2/3)  (132) 

<?ni  =  ai3sin(P)  +  a12cos(/3) 

in  which  /?  defines  the  principal  direction.  Matrix  failure  in  tension  is  given  by 

Onn 2  Ont2  Oni2  „  n  /i  oo\ 

+y  +  jr  >1  for  ann>0  (133) 

and  matrix  failure  in  compression  is  given  by 

Y~  +  yr  >1  for  onn<  0  (134) 
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4.2  Progressive  Material  Failure  Prediction 

Estimating  the  required  load  to  cause  the  next  material  failure  is  performed  by  calculating 
scale  factors  for  each  ply  in  each  element  to  precipitate  material  failure  using  selected  failure 
criteria  and  assuming  a  linear  load-deflection  response.  These  scale  factors,  a,  are  defined 
as  the  linear  multiplier  to  the  current  stresses  to  reach  the  critical  value  given  by 

oCT  -  aa  (135) 

such  that  the  failure  index,  FI,  is  equal  to  unity  or 

FI  =  =  1  (136) 

In  the  general  tensorial  criteria  described  previously,  this  leads  to  determining  roots  to  an 
nth  order  polynomial  of  the  form 


Snotn  +  Sn-\Oin  1  +  ...  +  So  —  0  (137) 

The  suite  of  failure  criteria  available  in  RESTRAN  are  limited  to  linear  and  quadratic  expres¬ 
sions  involving  stresses  and  material  strengths  such  that  the  highest  order  failure  criterion 
requires  only  the  roots  to  a  second-order  equation  of  the  form 

S20 ?  +  Sicx  -(-  Sq  —  0  (138) 

to  be  determined  to  yield  the  minimum  positive  multiplier,  a.  This  factor  is  calculated  for 
each  ply  in  each  element  and  is  stored  in  an  external  file.  The  lowest  factor  to  cause  the 
next  material  failure  is  selected  as 

a  =  mm(oT,  o?-, o£)  (139) 

where  i  is  the  element  number,  j  is  the  element  ply  number,  and  k  is  the  failure  mode. 
To  speed  convergence,  a  factor  s  is  applied  to  the  determined  load  multiplier  to  provide  a 
scaling  slightly  above  the  minimum  predicted.  This  factor  may  be  input  to  override  the 
default  value  of  1.01.  These  factors  multiply  the  vector  of  input  loads  yielding  the  next  set 
of  applied  loads  as 

F  =  saP  (140) 

In  the  case  when  nonlinear  material  properties  are  specified,  the  linear  or  quadratic  extrapo¬ 
lation  used  to  determine  a  will  introduce  an  additional  nonlinear  effect  and  cause  the  estimate 
for  the  next  failure  load  to  vary  during  the  iterative  solution  of  the  nonlinear  equilibrium 
equations. 

4.3  Material  Damage  Models 

Once  ply  failure  is  predicted,  various  damage  laws  have  been  incorporated  into  RESTRAN  to 
account  for  material  degradation.  The  simplest  damage  law  reduces  all  lamina  properties  to 
zero  regardless  of  damage  mode.  Thus,  for  the  ith  ply,  a  null  constitutive  relation  is  assumed 
as 

[C]<  =  [0]  .  (141) 
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This  assumption  is  regarded  as  extreme,  and  more  specialized  property  reductions  models 
have  been  advanced  for  the  representation  of  ply  damage  in  the  attempt  to  refine  failure 
predictions. 


Extending  the  model  developed  by  Reddy  and  Pandey  [56]  to  three-dimensional  for  com¬ 
pressive  or  tensile  matrix  failure,  the  in-plane  transverse  modulus,  E2,  G2z,  and  Poisson’s 
ratios  in  the  damaged  ply  is  set  to  zero  while  all  other  elastic  constants  are  unchanged. 
This  results  in  the  ply  stiffness  matrix  for  the  ith  ply  given  by 

'  Ei  0.0  0.0 

0.0  0.0  0.0 

r„1  _  0.0  0.0  £3 

1  Ji_  0.0  0.0  0.0 

0.0  0.0  0.0 

_  0.0  0.0  0.0 

Similarily,  for  fiber  failure,  E\,  G 12,  G13,  and  /Zy,  =  0,  yielding 

'  0.0  0.0  0.0 

0.0  E2  0.0 

rrl  _  0.0  0.0  Ez 

M<“  0.0  0.0  0.0 

0.0  0.0  0.0 

0.0  0.0  0.0 


A  more  specialized  damage  law  is  presented  by  Chang  and  Chang  [18].  This  law  is  based  on 
an  interactive  degradation  model  using  the  Tsai-Wu  failure  criterion  expressed  as 

FI  =  Fi<Ji  +  FijOiOj  (144) 

where  FI  is  the  failure  index.  Specific  failure  modes  are  determined  by  calculating  the  failure 
index  using  selected  stress  components.  The  implementation  in  RESTRAN  generalizes  this 
approach  to  any  selected  quadratic  failure  criterion.  A  stiffness  reduction  coefficient,  a,  is 
utilized  to  degrade  material  properties.  This  coefficient  is  experimentally  determined  and 
is  assumed  as  a  material  property  constant  for  the  composite  material  system  being  used. 
Fiber  breakage  is  assumed  if  max{FI(ai)}  is  due  to  o\.  The  damage  law  for  this  mode  is 
given  by 


Ei  =  aEi 

G\z  =  otG\z 

G\2  —  olGu  (145) 

P13  =  aviz 

P12  =  CM/12 
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with  all  other  moduli  left  unaltered.  Matrix  cracking  is  assumed  if  max{FI(oi)}  is  due  to 
<j2  or  oq,  the  corresponding  damage  law  for  this  mode  is  given  by 


damage  law  for  this  mode  is  given  by 
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For  partial  failure,  the  material  constitutive  matrix  will  generally  be  rank  deficient,  which 
will  lead  to  an  nth -order  singularity  during  inversion.  This  is  processed  through  an  element- 
level  condensation  of  the  indefinite  moduli  prior  to  the  inversion.  The  reduced  constitutive 
matrix  is  subsequently  expansed  to  the  original  material  matrix  dimension  and  is  then  used 
in  forming  element  stiffness  coefficients.  The  resulting  element  stiffness  matrix  will  possess 
zero-energy  modes  in  addition  to  those  associated  with  rigid  body  modes.  After  assembly 
of  elements  into  a  global  system,  any  degrees  of  freedom  associated  with  zero  stiffness  are 
statically  condensed  out  of  the  system  prior  to  analysis. 


The  failure  and  damage  models  built  into  RESTRAN  reflect  both  those  in  common  use 
and  those  that  are  the  result  of  current  research  into  modeling  failure  mechanisms  in  com¬ 
posites.  However,  RESTRAN  incorporates  a  standardized  subroutine  interface  described 
in  the  user’s  manual  [31]  which  allows  any  additional  material  failure/damage  model  to  be 
applied  through  the  creation  of  a  user-defined  subroutine  linked  into  the  RESTRAN  code. 


5  Structural  Failure  Modes 

Structural  members  exhibit  a  variety  of  buckling  responses  to  applied  loads.  Long,  slender 
members  such  as  beams  or  columns  tend  to  displace  in  a  simple  fundamental  mode  with  a 
single  wavelength  deflection  pattern.  Plate- like  components  demonstrate  various  numbers  of 
longitudinal  and  transverse  wavelength  deflection  patterns  which  depend  on  geometry  and 
material  aspect  ratios.  Thick  composite  laminates  with  delaminations  show  global,  local,  and 
mixed  instability  modes  as  depicted  in  Figure  11.  In  complex  geometries,  instability  modes 
can  preclude  any  clear  categorization.  The  algorithmic  treatment  of  delamination  buckling 
in  RESTRAN  allows  the  specification  of  any  number  of  delaminations  to  contribute  to.  local 
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and  mixed  instability  modes.  Each  delamination  may  be  defined  with  arbitrary  configura¬ 
tion  and  location  at  ply  interface  planes.  The  growth  of  delaminations  prior  to  buckling 
is  currently  not  calculated  in  RESTRAN.  This  exclusion  is  based  in  part  on  the  necessary 
limitation  of  the  current  developmental  effort,  together  with  the  observed  behavior  of  actual 
delminated  composites  undergoing  buckling.  Due  to  imperfections,  misalignments,  and  loss 
of  load  symmetry,  predictions  of  delamination  growth  prior  to  or  concurrent  with  buckling 
using  an  idealized  mathematical  model  tend  to  give  way  to  growth  occurring  subsequent  to 
buckling  and  being  a  phenomenon  in  the  post-buckled  regime  [11]. 

For  a  computationally  viable  analysis,  buckling  failure  is  predicted  through  a  linear  eige- 
nanalyis  to  obtain  critical  load  multipliers  to  cause  local  sublaminate  and  overall  instability. 
As  will  be  discussed  in  the  following  sections,  issues  arise  in  modeling  contact  constraint 
effects  and  in  representing  post-buckled  material  behavior.  To  fully  simulate  these  phenom¬ 
ena,  an  iterative  nonlinear  large  displacement  solution  procedure  is  required  which  would 
present  an  impractical  computational  cost  in  processing  multiple  delaminations.  Thus,  var¬ 
ious  simplifications  have  been  incorporated  to  provide  a  tractable  analysis. 

The  interpretation  of  buckling  modes  has  been  automated  for  delaminations  defined  along 
ply  interface  planes.  For  more  complicated  states  of  delamination  damage,  RESTRAN  incor¬ 
porates  user-defined  subroutine  interfaces  to  allow  problem-specific  interpretation  of  buckling 
failure  modes  for  complicated  delamination  profiles  in  three-dimensional  geometries. 


Figure  11.  Global,  local,  and  mixed  buckling  modes  in  thick  laminates. 
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5.1  Multiple  Delaminations 

The  analysis  of  multiple  delaminations  in  composites  has  not  been  widely  investigated  [60]. 
In  RESTRAN,  delamination  damage  is  idealized  as  an  array  of  planes  of  discontinuity;  this 
is  schematically  depicted  in  Figure  12.  In  low-velocity  impacts  of  plate-like  geometries,  the 
location  of  internal  delaminations  is  highly  dependent  on  the  relative  orientation  of  ply  pairs 
through  the  layup.  Other  effects  such  as  bending  and  stress  wave  interactions  tend  to  con¬ 
centrate  the  creation  of  fracture  surface  or  delamination  planes  towards  the  opposite  face  of 
the  panel  [61-66].  In  high  velocity,  through-penetration  type  impact  events,  delaminations 
may  be  uniformly  distributed  through  the  laminate  thickness  concentric  about  a  removed 
hollow  core. 


The  stability  analysis  of  multiply-delaminated  laminates  is  complicated  by  the  number  of 
potential  instability  modes.  If  delaminations  are  simply  modeled  as  planes  of  element  dis¬ 
continuity  and  a  linear  buckling  analysis  is  performed  with  out  enforcing  contact  constraints, 
predicted  buckling  modes  may  include  both  physical  modes  and  nonphysical  modes  in  which 
the  instability  failure  of  interior  sublayers  deform  into  surrounding  laminate  material.  This 
behavior  is  depicted  in  Figure  13. 


Figure  13.  Example  of  a  nonphysical  sublaminate  buckling  mode. 
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With  the  possible  presence  of  numerous  delaminations  with  arbitrary  location  and  configura¬ 
tion,  several  assumptions  have  been  used  in  the  development  of  a  general  analysis  capability 
in  RESTRAN  to  make  the  analysis  of  multiple  delaminations  tractable.  The  first  assumption 
is  that  each  sublaminate  buckling  instability  can  be  considered  individually.  Experimental 
results  have  shown  that  with  the  presence  of  multiple  delaminations,  instability  behavior  is 
typically  dependent  on  the  most  critical  delamination,  while  the  other  delaminations  remain 
closed  [11].  This  procedure  assumes  that  during  instability  of  a  particular  delamination,  all 
other  delamination  planes  have  infinite  frictional  contact  and  are,  as  a  consequence,  con¬ 
densed  out  of  the  model.  This  is  depicted  in  Figure  14.  This  process  is  repeated  sequentially 
for  all  delaminations  in  the  model,  and  the  delamination  with  the  lowest  critical  load  is 
assumed  to  exhibit  the  next  failure  due  to  buckling.  The  task  reduces  to  automating  the 
analysis  of  single  buckling  mode  shapes  and  ultimately  determining  which  particular  sublam¬ 
inate  layer  is  most  critical  to  instability  and  determining  which  elements  are  involved  in  the 
buckled  sublaminate.  This  procedure  eliminates  the  difficulty  of  automatically  screening  out 
nonphysical  buckling  modes  and  replaces  simultaneous  buckling  of  multiple  delaminations 
with  a  sequence  of  individual  local  instability  failures.  The  sequential  processing  of  individ¬ 
ual  local  delamination  failures  will  provide  overall  estimates  for  the  maximum  applied  load  to 
cause  local  failure  via  sublaminate  buckling.  In  a  study  of  local  buckling  in  laminated  com¬ 
posites,  the  presence  of  additional,  symmetrically  located  delaminations  have  only  a  slight 
effect  on  the  predicted  buckling  load  over  those  predicted  with  other  delaminations  removed 
from  the  model  [12].  A  schematic  of  sequential  buckling  failure  is  depicted  in  Figure  15. 


Figure  14.  Condensation  of  multiple  delaminations. 


Figure  15.  Sequential  buckling  failure. 

Although  sequentially  processing  individual  delaminations  removes  the  possibility  of  predict¬ 
ing  embedded  sublaminates  to  undergo  buckling  deformation  and  deforming  through  other 
layers,  as  depicted  in-  Figure  13,  contact  constraints  must  be  imposed  on  isolated  delami¬ 
nations  to  correct  for  local  interpenetration  of  the  buckling  mode,  as  shown  in  Figure  16. 
RESTRAN  incorporates  an  effective  iterative  procedure  for  satisfying  contact  or  compata- 
bility  constraints  on  the  buckling  mode  shape.  This  feature  is  explained  in  detail  in  the 
following  subsection. 


Figure  16.  Local  compatability  violation  in  layer  buckling. 


5.2  Algorithmic  Assessment  of  Local  Buckling  Failure 

Delaminations  are  assumed  to  exist  on  distinct  planes  which  may  be  of  arbitrary  shape  and 
extent.  Each  delamination  is  defined  by  a  declared  set  of  nodes  through  which  the  fracture 
surface  extends.  For  each  node  set,  RESTRAN  automatically  generates  a  set  of  coincident 
nodes  and  redefines  the  element  connectivity  to  provide  a  kinematic  freedom  of  motion  be¬ 
tween  the  upper  and  lower  surfaces.  Each  element  is  accessed,  and  the  nodes  forming  the 
element  faces  F2,  F3,  and  F6  (as  defined  in  Figure  17)  are  scanned  for  inclusion  in  the  delam¬ 
ination  node  set.  Once  determined,  the  elements  are  scanned,  and  those  containing  nodes  on 
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the  opposing  F4,  F\,  or  F$  faces,  which  are  a  subset  of  the  nodes  defining  the  delamination, 
are  assigned  the  coincident  nodes.  The  resulting  placement  of  coincident  nodes  is  depicted 
in  Figure  18. 


Figure  17.  Designation  of  element  faces. 


Ni  +  Incr  Nj  +  Incr  Nk  +  Incr  Ni  +  Incr 


Figure  18.  Generation  of  coincident  nodes. 

The  global  elastic  and  differential  stiffness  matrices  are  then  formed  containing  all  generated 
coincident  node  planes  simulating  the  required  delaminations.  A  sequential  processing  of 
each  delamination  is  performed,  and  a  condensation  is  performed  to  temporarily  remove  all 
but  the  current  delamination.  An  eigenanalysis  of  the  condensed  global  system  is  then  in¬ 
voked.  Each  fundamental  mode  shape  is  subjected  to  an  automated  interpretation.  The  first 
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test  determines  the  type  of  buckling  mode.  Simple  global  buckling  is  identified  if  no  delam¬ 
inations  have  been  input  or  after  all  delaminations  have  been  failed.  A  mixed  global-local 
buckling  mode  is  detected  if  the  maximum  normalized  modal  displacement  occurs  outside 
the  domain  of  the  delamination  plane.  This  mode  is  predicted  when  the  ratio  between  the 
maximum  modal  deflection  in  the  delamination  plane  and  the  maximum  overall  deflection 
is  less  than  an  input  tolerance  or 


«/■&  <  Tol  (148) 

If  these  modes  are  not  exhibited,  the  default  mode  is  local  sublaminate  buckling.  These 
buckling  modes  are  depicted  in  Figure  11. 

Global  or  mixed  global-local  buckling  modes  are  assumed  to  cause  catastrophic  failure  if 
the  corresponding  critical  load  is  less  than  the  predicted  load  to  cause  material  failure  at  the 
current  progressive  failure  cycle. 

For  local  sublaminate  buckling  modes,  an  iterative  procedure  can  be  invoked  to  satisfy 
contact  constraints.  The  automated  analysis  of  the  buckling  mode  is  based  entirely  on  the 
mode  shape  which  is  arbitrary  in  magnitude  and  sign.  An  estimate  for  the  amplitude  of 
the  buckling  mode  within  a  linear  solution  scheme  can  be  made  by  assuming  the  area  of  the 
delaminated  section  remains  constant  and  that  the  membrane  stress  in  the  buckled  laminate 
is  the  same  as  the  buckling  stress  [6].  For  assessing  contact,  however,  the  actual  magnitude 
is  not  needed,  and  the  pattern  of  the  normalized  eigenmode  is  used  to  determine  modal 
displacements  and  modal  stresses  and  strains.  Iterations  begin  by  first  checking  whether  the 
buckling  mode  is  demonstrating  a  physical  opening  mode  or  nonphysical  closing  mode,  as 
shown  in  Figure  19. 


OPENING  MODE 

Figure  19.  Opening  and  closing  buckling  modes. 
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If  a  closing  mode  is  detected,  the  arbitrary  signs  of  the  modal  displacements  are  flipped.  A 
test  is  performed  to  assess  the  degree  of  interpenetration  of  the  buckling  layer  into  the  upper 
or  lower  substrate.  The  inital  processing  of  interpenetrating  delaminations  is  to  remove  the 
coincident  nodes  where  contact  conditions  are  violated.  This  has,  as  a  first  approximation, 
the  effect  of  fitting  an  opening  half  wavelength  to  the  delamination  plane.  A  minimum  energy 
surface  may,  however,  extend  beyond  this  border,  and  further  computations  are  performed 
to  refine  the  inital  estimate  of  constrained  layer  buckling.  This  is  performed  by  repeating 
the  eigenanalysis  with  the  augmented  connectivity  and  utilizing  the  recovered  mode  shape 
which,  although  arbitrary  in  magnitude  and  sign,  provides  the  necessary  qualitative  infor¬ 
mation.  After  testing  for  interpenetrating  nodes  and  accumulating  these  in  an  exclusion  set, 
prior  nodes  included  in  this  set  are  tested  for  possible  release  due  to  predicted  normal  tensile 
strains  at  the  node  location.  This  estimate  is  based  on  tensile  stresses  associated  with  nodes, 
or  equivalently,  nodal  displacements  above  and  below  the  delamination  plane  moving  apart, 
thereby  generating  tensile  strains  normal  to  the  plane  which  act  to  open  that  portion  of  the 
delamination.  These  nodes  are  restored  in  the  delamination  node  set  and  the  linear  analysis 
rerun.  The  nodes  used  for  estimating  modal  strains  about  a  particular  coincident  node  are 
shown  in  Figure  20. 


Figure  20.  Nodes  used  to  calculate  normal  modal  strain. 


The  coordinates  Xi  and  modal  deflections  6i  are  then  used  to  calculate  a  measure  of  normal 
strain  as  a  stretch  factor  given  by 


£77.  — 


Si  -  s0  [zU(xt  +  -  sr)2f2  -  [zh(xt  -  *r)2]1/2 
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(149) 


If  the  criteria  en  >  Tol  is  satisfied,  the  coincident  node  previously  removed  through  con¬ 
densation  is  reintroduced.  Other  approaches,  such  as  using  iterative  conjugate  gradient 
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methods,  have  been  applied  to  satisfy  buckling  layer  contact  with  linear  constraints  [67,68]. 
The  present  approach  within  the  context  of  linear  bifurcation  buckling  provides  a  similar 
solution  in  a  small  number  of  iterations.  To  perform  an  iterative  contact  solution  for  ac¬ 
ceptable  mode  shapes,  user  input  to  RESTRAN  specifies  the  maximum  number  of  iterations 
together  with  a  tolerance  which  is  a  normalized  measure  of  overall  contact  above  which  the 
compatibility  is  acceptable.  This  measure  is  given  by 

1  -  ^  >  T„(  (150) 

where  Njp  is  the  number  of  interpenetrating  nodes,  and  Np  is  the  total  number  of  nodes  in 
the  delamination  set.  If  convergence  is  not  obtained,  the  mode  is  assumed  impossible  and  is 
excluded  from  consideration  in  the  current  analysis  cycle. 

For  a  physically  acceptable  local  sublaminate  buckling  mode,  the  next  procedure  is  to  de¬ 
termine  whether  the  sublaminate  layer  undergoing  buckling  is  along  the  positive  or  negative 
normal  to  the  delamination  plane.  As  shown  in  Figure  21,  the  nodes  defining  the  element 
face  and  contained  in  the  delamination  node  set  are  used  to  form  two  local  vectors,  Vi  and 
V2,  which  define  the  local  positive  normal  to  the  delamination  surface  as 

Vi  x  V2 
n  \Vi  x  V2| 


Figure  21.  Definition  of  positive  normal  to  delamination  plane. 

The  buckling  surface  is  determined  by  forming  displacement  vectors  from  nodes  in  the  de¬ 
lamination  plane  to  the  upper  and  lower  surfaces  and  computing  vector  norms.  The  surface 
which  demonstrates  the  greatest  magnitude  of  displacement  then  determines  which  sublam¬ 
inate  is  undergoing  buckling.  Figure  22  shows  a  depiction  of  an  opening  delamination  with 
local  displacement  vectors  at  node  points. 
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Figure  22.  Opening  mode  in  delamination  buckling. 

The  direction  of  opening  with  respect  to  the  delamination  surface  normal  is  determined  by 
computing  the  angle  between  the  nodal  displacement  vectors  and  the  surface  normal.  The 
individual  angles  are  summed,  and  an  average  opening  angle  a  is  computed. 


a  = 


1  N 

—  $2 cos 
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n 

Wi  ' 

>i 
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The  direction  d  is  thus  determined  as 


(152) 


d  = 


+n 
— n 


Of  <  7t/2 
a  >  7t/2 


(153) 


In  the  event  of  a  centrally  located  delamination  (Figure  23)  wherein  both  surfaces  are 
moving  symmetrically  apart,  a  may  be  identically  equal  to  7r/2,  and  the  positive  normal 
direction  is  taken  by  default.  This  will  cause  the  proper  critical  applied  load  to  be  ac¬ 
counted,  but  only  one  of  the  layers  will  be  failed.  The  other  layer  will  remain  and  be  as¬ 
sessed  during  the  next  analysis  cycle.  Again,  this  suggests  that  sequencing  failures  through 
sublaminate  buckling  compared  with  attempting  to  account  for  simultaneous  failures  will 
alter  the  progression  of  failure  while  maintaining  an  accurate  estimate  for  failure  loads. 
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Figure  23.  Buckling  failure  in  a  centrally  located  delamination. 


As  each  delamination  is  processed,  a  continuous  update  is  made  keeping  account  of  the 
delamination  with  the  lowest  predicted  critical  buckling  load,  associated  mode  shape,  and 
direction  of  failure.  If  no  other  material  failures  are  predicted  at  a  lower  load  level,  the 
selected  critical  delamination  is  then  designated  as  the  next  failure  event.  Processing  the 
failed  delamination  consists  of  performing  a  search  through  the  sublaminate  undergoing 
buckling  and  determining  what  elements  are  involved  in  the  instability  failure.  Once  de¬ 
termined,  a  selected  damage  law  is  applied  to  each  of  the  elements.  An  additional  search 
is  performed  to  determine  if  any  other  delaminations  are  contained  in  the  buckled  sub¬ 
laminate.  This  is  depicted  in  Figure  24.  All  delaminations  involved  in  the  current  buck¬ 
ling  failure  are  then  removed  from  the  node  set  and  removed  from  the  model  by  elimi¬ 
nating  the  involved  coincident  nodes.  An  internal  node  renumbering  is  then  performed. 


Figure  24.  Failure  assessment  of  buckled  sublaminate. 
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5.3  Distributed  Delaminations  and  Simultaneous  Local  Buckling 

The  algorithms  incorporated  into  RESTRAN  for  automatically  processing  delamination 
buckling  can  be  used  to  analyze  more  complicated  delamination  configurations.  For  ex¬ 
ample,  a  traversing-type  delamination  which  progresses  along  several  ply  interfaces  can  be 
modeled.  As  shown  in  Figure  25,  the  user  must  initially  create  the  model  with  connectivity 
removed  between  elements  along  the  perpendicular  segments  where  the  delamination  jumps 
between  different  ply  interfaces.  The  remaining  portions  of  the  delamination  can  be  defined 
by  specifying  a  delamination  node  set.  Additionally,  as  shown  in  Figure  26,  delaminations 
do  not  have  to  be  similarily  oriented.  Each  independent  delamination  may  be  arbitrarily  po¬ 
sitioned  in  different  areas  of  the  structural  model.  The  only  restriction  is  that  each  separate 
delamination  must  form  a  planar  surface  for  proper  interpretation  of  buckling  motion.  To 
allow  the  capability  of  analyzing  general,  nonplanar  delamination  surfaces  or  the  simultane¬ 
ous  buckling  of  multiple  delaminations,  RESTRAN  provides  user-defined  subroutine  options 
to  interpret  buckling  in  non-plate-type  geometries.  Delaminations  are  defined  in  RESTRAN 
by  inputting  node  sets  which  are  used  to  automatically  generate  coincident  nodes  to  redefine 
element  connectivity.  Thus,  any  discontinuous  inner  surface  may  be  simulated,  and  multiple 
delaminations  may  be  considered  simultaneously  by  including  the  participating  nodes  in  the 
same  delamination  node  set.  A  user-defined  subroutine  which  is  compiled  and  linked  into 
the  RESTRAN  executable  must  then  be  provided  to  interpret  the  resulting  mode  shapes 
and  assess  associated  element  failures. 


Figure  25.  Modeling  traversing  delaminations. 
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5.4  Post-Buckled  Failure  Modeling 

The  stability  analysis  incorporated  into  RESTRAN  is  linear.  However,  because  local  buckling 
is  predicted  during  the  progressive  failure  analysis,  some  approximation  to  post-buckled  layer 
response  must  be  applied  while  failure  of  the  remaining  laminate  is  being  calculated.  Towards 
this  aim,  after  the  prediction  of  local  buckling  is  made,  the  coincident  nodes  defining  the 
failed  delamination  are  removed  through  condensation  from  the  model,  and  assumptions  are 
made  regarding  the  remaining  strength  of  the  elements  participating  in  the  buckling  mode. 
As  depicted  in  Figure  27,  different  structures  exhibit  a  range  of  post- buckled  load-deflection 
response.  As  shown,  a  plate-type  structure  can  exhibit  significant  further  load  carrying 
capability  with  possible  transition  to  higher-order  mode  shapes  with  an  effective  reduction 
in  material  stiffness.  To  approximate  post-buckled  material  behavior,  RESTRAN  permits  an 
inelastic  reduction  in  moduli  to  be  applied.  This  implies  a  permanent  reduction  in  properties; 
no  elastic  recovery  of  properties  is  assumed  possible  with  unloading  of  the  buckled  layer. 
Several  assumptions  on  post-buckled  material  behavior  may  be  specified  in  RESTRAN. 
In  the  extreme  case,  moduli  may  be  set  to  zero  to  remove  any  load  carrying  capability 
in  the  buckled  layer  and  thus  redistribute  the  applied  loads  to  the  remaining  structure. 
An  alternative  assumption  is  that  the  buckled  layer  can  further  absorb  load  -  and  exhibit 
possible  subsequent  or  additional  material  failure  modes  -  but  at  a  reduced  modulus.  Finally, 
the  intermediate  assumption  of  carrying  the  constant  load  in  the  post-buckled  sublaminate 
corresponding  to  the  critical  load  is  approximated  in  RESTRAN  through  the  option  of  a 
scaling  procedure  of  element  stiffness  coefficients.  This  scaling  for  maintaining  a  constant 
load  in  a  buckled  sublaminate  is  accurate  to  the  degree  that  the  resulting  layer  deformations 
at  subsequent  loading  remain  similar  to  the  pattern  existing  at  the  buckling  failure  load. 
This  may  be  formulated  by  considering  the  initial  element-level  equilibrium  expressed  as 

{F„}  =  [K]{U„)  (154) 
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where  {F0}  are  internal  forces  due  to  the  initial  set  of  external  applied  loads.  For  a  linear 
structure,  the  equilibrium  relation  at  a  critical  multiple  of  the  initial  loads  can  be  given  by 

aCr{F0}  =  [K]acr{U0}  (155) 

If  the  displacement  field  remains  similar  such  that  at  a  higher  load  scale  factor  a*  yielding  a 
set  of  displacement  {U;},  it  remains  true  that 

{U0}  «  -{UJ  (156) 

Oii 

Then,  for  any  applied  load,  a  scaling  of  the  element  stiffness  matrices  for  those  elements 
involved  in  the  buckled  sublaminate  given  by 


*cr{F0}  =  /?[K]{UJ  (157) 

will  fix  the  equilibrium  element  forces  at  the  level  at  which  buckling  was  predicted  with  j3 
given  by 

0  =  —  (158) 

Oii 

Or,  to  allow  increasing  loads,  /3  is  implemented  in  RESTRAN  as 

^(?-1)Ci+1  (159) 

where  Ci  may  be  varied  between  0  and  1.  Thus,  representing  post-buckled  behavior  in  a 
linear  analysis  is  highly  approximate,  and  the  selection  of  a  plausible  post-buckled  elastic 
response  is  important  to  guarantee  a  conservative  estimate  of  the  overall  residual  strength. 


Figure  27.  Post-buckling  load  carrying  capability. 
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6  Solution  Algorithms 

Four  different  analysis  procedures  have  been  incorporated  into  RESTRAN.  As  depicted 
schematically  in  Figure  28,  the  available  options  are  Prespass,  Linear  Static  Analysis,  Linear 
Buckling  Analysis,  and  Residual  Strength  Analysis. 


Figure  28.  RESTRAN  solution  control  options. 


6.1  Prepass 

Selecting  prepass  invokes  assembling  and  checking  the  inputted  model  after  which  the  pro¬ 
gram  is  terminated.  This  causes  numerous  tests  on  the  finite  element  model  to  be  performed 
which  are  written  to  the  standard  output  file.  The  internal  tests  performed  include  checks 
on  element  connectivity  and  geometry,  boundary  conditions,  applied  point  or  pressure  loads, 
and  material  property  assignments.  In  addition,  RESTRAN  provides  an  optional  feature 
which  creates  an  output  file  containing  graphics  information  to  view  the  input  model.  RE¬ 
STRAN  supports  output  formats  for  MATHEMATICA  [69]  or  TECPLOT  [70]  and  has  the 
option  for  a  user-defined  subroutine  to  format  graphics  data  as  required.  This  feature  is 
used  as  an  additional  test  to  check  for  possible  modeling  errors  prior  to  performing  a  de¬ 
tailed  residual  strength  analysis.  The  graphical  depiction  is  particularly  useful  in  complex 
geometries  which  may  have  been  created  using  general  preprocessing  codes  such  as  PATRAN 
[71],  or  using  the  node  and  element  generation  option  available  in  RESTRAN.  The  flowchart 
of  execution  in  Prepass  is  depicted  in  Figure  29. 


49 


Figure  29.  Prepass  flow  chart. 


6.2  Static  and  Buckling  Analysis 

The  single-pass  displacement  and  buckling  solutions  are  made  available  to  obtain  the  linear 
elastic  response  of  the  model  under  applied  loads.  Schematic  flow  charts  of  these  solution 
sequences  are  presented  in  Figures  30  and  31.  The  displacement  solution  yields  the  basic 
deformation  pattern  together  with  optional  stress  and  strain  output  at  ply  layers  or  Gauss 
points.  In  addition,  requesting  a  ply  failure  output  will  cause  additional  computations  of 
failure  indices  within  elements.  This  information  is  used  to  output  the  minimum  predicted 
scale  factor  to  the  applied  loads  to  cause  the  next  material  failure.  Within  the  context  of 
a  linear  static  analysis,  this  scale  factor  may  be  considered  as  a  residual  strength  measure 
based  on  an  overall  first-ply  failure  prediction.  The  buckling  solution  may  be  used  as  a  pre¬ 
liminary  check  on  single-layer  buckling  response  or  the  simultaneous  instability  of  multiple 
delaminations  as  a  comparison  to  running  a  residual  strength  analysis  in  which  multiple 
delaminations  are  automatically  analyzed  individually.  In  both  displacement  and  buckling 
solutions,  graphical  output  may  be  requested  that  is  generated  in  the  form  of  an  external  file 
which  may  be  directly  inputted  into  a  selected  graphics  program  such  as  MATHEMATICA 
or  TECPLOT  to  view  deformed  geometry  or  buckling  mode  shape. 
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ASSEMBLE  GLOBAL  STIFFNESS 
MATRIX  AND  FORCE  VECTOR 


NEL  NEL T  _1 

[K]c=BK4=IJJJ[G.]  [H,]  [G,] 

i  =  l  i  =  1 

NODE  FACE 
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i  =  1  i  =  1 


Figure  30.  Flow  chart  of  static  analysis. 
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Figure  31.  Flow  chart  of  buckling  analysis. 
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6.3  Residual  Strength  Analysis  Algorithm 

In  representing  realistic  structural  geometries  and  support  and  load  conditions,  a  viable  al¬ 
gorithmic  treatment  in  simulating  material  behavior  and  failure  modes  under  increasing  load 
is  critical  to  the  success  of  the  analysis.  In  real  structures,  the  progression  of  local  failure 
modes  are  not  independent,  but  tend  to  be  coupled  and  occur  simultaneously.  A  sequen¬ 
tial  algorithm,  however,  to  be  tractable,  requires  that  these  effects  be  considered  separately. 
The  important  structural  behaviors  simulated  by  RESTRAN  are  nonlinear  material  moduli, 
material  ply-level  failure  differentiated  into  various  modes,  local  instability  of  delaminated 
sublaminates  and  global  buckling  of  the  entire  model,  and  approximate  post-buckled  mate¬ 
rial  response.  The  primary  variable  is  the  scale  factor  applied  to  the  vector  of  initial  loads 
to  determine  a  sequence  of  failure  in  a  combined  incremental  and  iterative  fashion.  The  esti¬ 
mated  scale  factor  for  causing  material  failure  is  obtained  from  linear  extrapolations  for  each 
ply  using  selected  failure  criteria.  The  scale  factor  for  minimum  buckling  loads  are  obtained 
from  eigenanalysis.  The  effects  of  material  nonlinearity  due  to  stress/strain  relationship 
and  post-buckled  behavior  of  delaminated  layers  are  considered  separately  from  failure  pre¬ 
diction.  Thus,  at  the  ith  failure  cycle,  iterations  are  first  performed  to  converge  both  the 
incremental  displacements  and  the  minimum  positive  buckling  eigenvalue.  At  convergence, 
the  load  multipliers  are  compared,  and  the  minimum  dictates  whether  the  next  failure  event 
is  predicted  to  be  a  material  mode  or  a  buckling  mode.  For  buckling  failure,  the  element 
set  involved  is  computed,  selected  post-buckled  failure  assumptions  imposed,  and  the  failed 
delamination  is  condensed  out  of  the  model.  For  material  failure,  iterations  are  performed  at 
fixed  load  to  update  element  properties  due  to  applied  damage  laws,  stress  redistribution  is 
calculated,  and  additional  element  failures  are  processed.  Convergence  is  established  when 
no  additional  element  failures  are  predicted  at  the  current  load.  A  flowchart  of  the  residual 
strength  prediction  algorithm  is  contained  in  Figure  32.  The  residual  strength  solution  anal¬ 
ysis  combines  a  number  of  available  options.  These  options  modify  the  solution  algorithm  to 
assess  assumed  modes  of  failure  or  to  tailor  the  analysis  for  specialized  failure/ damage  laws 
or  complicated  geometries.  Table  3  shows  the  analysis  options. 


Table  3.  RESTRAN  solution  options. 


•  Predict  progressive  material  failures  only. 

•  Predict  sequential  buckling  failures  only. 

•  Predict  combined  material  and  instability  failures. 

•  Apply  standard  material  failure  criteria. 

•  Apply  user-defined  material  failure  criteria. 

•  Apply  standard  material  damage  laws. 

•  Apply  user-defined  material  damage  laws. 

•  Access  user-defined  subroutine  to  interpret  complex  buckling 
modes  such  as  simultaneous  multiple  delamination  buckling. 

•  Process  nonlinear  elastic  material  behavior. 

•  Apply  specific  post-buckled  material  behavior. _ 
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Figure  32.  Flow  chart  of  residual  strength  solution  algorithm. 
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7  Ultimate  Failure  Prediction 


The  residual  strength  of  a  damaged  laminated  composite  is  obtained  when  catastrophic  or 
ultimate  failure  is  predicted.  During  execution  of  the  solution  sequence  in  which  incremental 
failures  are  assessed,  a  running  account  of  the  maximum  scale  factor  to  the  applied  loads  is 
maintained.  A  typical  failure  sequence  may  initially  demonstrate  low  loads  to  cause  local 
failures  such  as  near-surface  delamination  buckling  or  material  failure  in  an  area  of  high 
stress  concentration.  Subsequent  failures  may  involve  simultaneous  buckling  of  thick  sub¬ 
laminate  ply  groups  or  large  regions  of  material  damage.  Final  failure  may  be  precipitated 
through  global  buckling  of  remaining  layers  or  through  a  catastrophic  cascade  of  element 
material  failures  as  loads  are  redistributed  at  each  iteration.  Specific  checks  performed  in 
RESTRAN  are  listed  in  Table  4. 

Table  4.  Tests  for  catastrophic  failure. 


•  Existence  of  ’soft’  deformation  or  rigid  body  displacement  modes. 

•  Prediction  of  global  or  mixed-mode  buckling. 

•  Stiffness  loss  exceeding  90%  due  to  material  and  structural  failure. 

•  Failure  at  nodes  at  which  external  loads  have  been  applied. 


If  the  lowest  load  multiplier  in  a  particular  analysis  cycle  predicts  global  buckling,  no  fur¬ 
ther  load-carrying  capability  in  the  post-buckled  regime  is  assumed,  and  ultimate  failure 
is  predicted.  As  element  material  properties  are  degraded  due  to  material  or  local  buckling 
failure,  possible  fragmentation  is  assessed  by  testing  for  positive  indefiniteness  by  performing 
an  LrDL  decomposition.  Related  to  this,  if  the  percentage  of  individual  degrees  of  freedom 
with  no  associated  stiffness  exceeds  90%,  ultimate  failure  is  assumed,  regardless  if  the  re¬ 
maining  model  is  unfragmented.  Finally,  if  failure  has  occurred  at  nodal  degrees  of  freedom 
to  which  external  loads  are  applied,  total  failure  is  assumed  because  subsequent  scalar  load 
multipliers  lose  definition  as  the  initial  vector  of  applied  loads  has  been  changed.  An  input 
parameter  is  available  to  preclude  elements  associated  with  applied  loads  from  exhibiting 
failure. 


8  Computer  Implementation  of  RESTRAN 

RESTRAN  is  written  in  FORTRAN  77.  In  the  development  of  the  various  algorithms,  no 
special  features  dependent  on  specific  computer  platforms  have  been  exploited  to  speed  exe¬ 
cution  in  order  to  guarantee  portability  of  the  code.  The  finite  element  basis  of  RESTRAN 
naturally  leads  to  the  generation  of  large,  sparse,  banded  matrices  in  representing  the  global 
elastic  and  differential  stiffness  properties  of  the  complete  model.  Internal  algorithms  have 
been  created  to  process  matrix  storage  modes  in  various  formats.  These  include  full  matrix 
storage  for  small  problems,  half-bandwidth  storage  format  (default),  and  out-of-core  storage 
for  large  problems  that  exceed  internal  RAM  memory  capacity  or  user  stacksize  limits  and 
require  most  data  storage  to  be  placed  in  external  files.  These  features  were  incorporated  as- 
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suming  that  not  all  platforms  on  which  RESTRAN  would  be  installed  would  support  internal 
virtual  memory  swap  operations.  Setting  the  internal  memory  parameters  and  selecting  the 
algorithmic  path  are  discussed  in  the  RESTRAN  User’s  Manual  [31].  The  choice  of  internal 
operation  is  entirely  dependent  on  execution  speed,  which  is  determined  by  the  amount  of 
data  that  can  be  held  in  high-speed  in-core  memory.  Thus,  the  out-of-core  solution  mode  is 
slowest  due  to  the  high  I/O  overhead  of  transfering  data  between  core  and  external  files.  Al¬ 
though  different  algorithmic  approaches  are  sometimes  used  for  the  different  memory  storage 
modes  to  perform  various  operations,  excluding  I/O  operations,  they  are  computationally 
competitive  in  terms  of  operation  counts. 

After  parsing  the  input  file,  element  connectivity  information  is  used  in  a  modified  Sloan’s 
bandwidth  minimization  procedure  [72]  to  obtain  internal  equation  numbers.  The  band¬ 
width  is  determined  as  the  maximum  of  the  bandwidth  resulting  from  considering  only  the 
initial  input  connectivity  and  from  that  determined  by  considering  the  coincident  nodes  due 
to  all  input  delamination  planes.  These  operations  are  performed  in  a  preface  section  of 
the  RESTRAN  code  which  liberally  allocates  memory  to  various  arrays  used  in  input  file 
processing.  After  basic  model  characteristics  are  determined,  such  as  the  number  of  degrees 
of  freedom,  nodes,  elements,  and  connectivity  as  represented  by  the  bandwidth,  the  internal 
memory  is  reconfigured  such  that  array  dimensions  are  streamlined  to  the  required  job  size. 

The  global  elastic  and  differential  stiffness  matrices  are  formed  using  a  direct  assembly 
approach.  For  in-core  operation,  the  entire  global  stiffness  matrix  is  accumulated  before 
being  written  out  to  disk.  In  an  out-of-core  operational  mode,  blocks  of  eight  rows  of  the 
global  stiffness  matrix  are  read  in  and  element  contributions  are  added  based  on  the  internal 
equation  numbering  of  element  degrees  of  freedom  (DOF).  This  block  is  then  written  to  an 
external  file  and  the  process  repeated  until  all  stiffnesses  have  been  assembled. 

During  execution,  degrees  of  freedom  are  removed  from  the  active  analysis  set  through 
applied  zero  displacement  boundary  conditions,  failed  DOF’s  at  which  material  failure  has 
reduced  the  stiffness  to  zero,  and  eliminated  DOF’s  associated  with  coincident  nodes  removed 
when  failed  delaminations  are  removed  from  the  model.  The  condensation  procedure  utilizes 
a  direct  row/column  elimination  on  the  global  matrices.  Disk  space  utilization  requires  that 
in  addition  to  the  elastic  and  differential  stiffness  terms  stored  externally,  decompositon  and 
reduction  of  these  data  blocks  are  written  to  separate  files.  This  leads  to  the  requirement 
that  four  global  matrices  need  to  be  stored  externally  in  scratch  files;  thus,  care  must  be 
given  to  configuring  an  adequate  size  of  the  disk  partition  in  which  RESTRAN  is  to  be 
executed.  Internally,  as  soon  as  a  matrix  is  no  longer  needed,  the  external  file  in  which  it  is 
resident  is  closed  to  free  up  disk  space. 

The  solution  for  displacements  is  obtained  using  different  equation  solvers  depending  on 
storage  mode.  For  full  matrix  storage  a  Gaussian  elimation  scheme  with  full  pivoting  is  em¬ 
ployed.  For  band  storage,  the  sparse  symmetric  banded  solver  routines  DPBCO  and  DPBSL, 
available  from  the  LINPACK  [73]  linear  algebra  library,  are  used.  The  solution  is  based  on 
a  Cholesky  decomposition  followed  by  a  back  substitution  using  the  triangular  Cholesky 
factor.  A  variation  of  this  algorithm  is  used  for  out-of-core  solution  and  is  implemented 
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in  RESTRAN  based  on  the  program  SESOL  presented  by  Wilson  et  al.  [74].  This  solver 
calculates  the  maximum  equation  block  size  that  can  be  held  in  core.  Upon  completing  the 
decomposition  of  each  block,  the  solution  vector  is  formed  by  sequentially  accessing  each 
block  contained  on  direct  storage.  The  minimum  memory  requirements  of  this  routine  is  for 
at  least  two  equations  being  held  in  memory  per  equation  block.  The  number  of  blocks  is 
dictated  by  the  amount  of  high-speed  RAM  memory  available. 

Eigenanalysis  using  full  and  band  storage  modes  can  be  performed  using  a  subspace  it¬ 
eration  method  detailed  in  section  2.4.  The  current  implementation  uses  the  implicitly 
restarted  Arnoldi  iteration  method  for  band  and  out-of-core  storage  modes.  A  robust  im¬ 
plementation  of  this  method  is  contained  in  the  routines  DSAUPD  and  DSEUPD  available 
in  the  ARPACK  distribution  [75].  In  turn,  these  routines  are  dependent  on  basic  operations 
and  algorithms  provided  by  the  BLAS  and  LAPACK  libraries  [76,77].  As  implemented  in 
RESTRAN,  minimum  memory  requirements  are  approximately  12  x  Nj,  where  iVj  is  the 
order  of  the  eigenproblem.  In  forming  iteration  vectors,  (j>,  it  is  required  to  repeatedly  solve 
a  linear  system  of  the  form 


{</>}  =  [K  —  A[Kq]]_1{KV}  (160) 

where  V  is  an  iteration  vector.  This  is  accomplished  using  an  appropriate  solver  discussed 
above  and  constitutes  the  major  computational  expense  in  performing  the  eigenanalysis. 

A  single  large  working  array  is  used  in  RESTRAN  to  perform  computations.  Total  memory 
requirements  for  in-core  storage  modes  are  essentially  dictated  by  the  size  of  this  array.  In 
a  full  matrix  storage  mode  used  for  small  problems,  core  memory  is  essentially  ND  x  ND 
double  precision  words,  where  ND  is  greater  or  equal  to  the  number  of  input  DOF.  For  band 
storage,  in-core  memory  is  on  the  order  of  No  x  Mb,  where  Mb  is  the  half  bandwidth  of 
the  global  matrix.  Finally,  using  the  out-of-core  algorithm  reduces  the  size  of  the  working 
array  such  that  the  collective  size  of  the  numerous  other  small  arrays  used  in  RESTRAN 
become  dominant  in  setting  the  program  memory  size.  In  the  current  program  release,  this 
requirement  is  approximately  equal  to  21  x  Njy- 


9  Numerical  Studies 

An  initial  set  of  benchmark  problems  were  solved  to  assess  the  accuracy  and  convergence 
properties  of  the  hexahedral  element  incorporated  into  RESTRAN,  and  to  suggest  a  degree 
of  discretization  required  to  adequately  simulate  local  laminate  behavior.  Of  particular  in¬ 
terest  is  in  quantifying  the  effect  of  the  incompatible  modes  incorporated  into  the  Pian-Tong 
element  in  accurately  supporting  bending  behavior  within  a  three-dimensional  continuum 
element  formulation.  The  Pian-Tong  element  is  well  known  to  be  insensitive  to  geomet¬ 
ric  distortion  and  has  a  superior  accuracy  over  displacement-based  elements  in  predicting 
stresses.  These  features  yield  an  ideal  element  for  modeling  three-dimensional  stress  states 
for  local  failure  prediction,  and  thereby  yield  accurate  differential  stiffnesses  which  are  evi¬ 
denced  by  the  rapid  convergence  of  buckling  load  predictions.  A  first  illustration  of  the  basic 
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element  performance  is  the  prediction  of  Euler  column  buckling  loads.  Three  models  shown 
in  Figure  33  were  selected  using  3,  5,  and  10  elements,  repectively.  The  column  dimensions 
are  10  x  1  x  1,  and  material  properties  are  given  by  E  =  1.0E9  and  u  =  0.3.  The  results 
presented  in  Table  5  show  a  rapid  convergence  to  the  exact  solution. 


Figure  33.  Models  for  Euler  column  buckling  load  determination. 


Table  5.  Convergence  of  column  buckling  loads. 


Model 

A 

%  error 

3-Element 

2187.00 

6.363 

5-Element 

2099.67 

2.126 

10-Element 

2060.65 

0.218 

Exact 

2056.17 

The  modeling  of  plates  using  three-dimensional  solid  elements  has  been  demonstrated  to 
generate  accurate  results  with  large  element  aspect  ratios  while  representing  full  three- 
dimensional  elasticity  [78-81].  To  illustrate  the  capability  of  the  current  three-dimensional 
continuum  element  to  predict  plate-like  behavior,  a  simply  supported  square  plate  was  se¬ 
lected  with  dimensions  10  x  10  x  0.1  and  isotropic  material  properties  given  by 

E  =  1.0000£9 

G  =  0.3846JE'9  (161) 

v  =  0.3 
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This  plate  was  analyzed  statically  to  determine  maximum  deflections  due  to  a  uniform  ap¬ 
plied  surface  pressure,  and  instability  was  analyzed  by  considering  the  plate  subjected  to 
end-loading  in  compression  and  determining  the  fundamental  buckling  load.  The  geometry 
and  loading  conditions  are  depicted  in  Figure  34.  The  convergence  of  the  solution  with  in¬ 
creasing  discretization  is  presented  in  Table  6. 


Figure  34.  Geometry  and  loading  of  a  flat  plate. 

Table  6.  Convergence  study  for  an  isotropic  plate  model. 


Model* 

Wmax  i 

A 

2  x  2f° 

-.0003707 

5.0569E4 

4  x  4f“ 

-.0004391 

3.8996E4 

4  x  4sfc 

-.0004431 

3.6791E4 

8  x  8s6 

-.0004441 

3.6252E4 

Exact 

-.0004443 

3.6152E4 

a  A  full  plate  was  used  in  generating  results. 

6  A  symmetric  quarter  plate  was  used  in  generating  results. 


To  assess  the  effect  of  aspect  ratio  in  the  solution  for  plate  deflections,  the  isotropic  plate 
used  above  was  analyzed  with  different  thickness  ratios.  The  solutions  are  presented  in  Ta¬ 
ble  7  and  are  compared  with  the  exact  solution  for  a  thin  plate  in  plane  stress  given  by 
Timoshenko  and  Woinowsy-Krieger  [82]  by 

Wmax  =  4.433  x  10"4  (J^)  (162) 
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Table  7.  Effect  of  element  aspect  ratio. 


a/h 

Wmclx/WcPT 

10 

1.1293 

100 

1.0017 

1000 

1.0000 

10000 

0.9194 

100000 

0.0035 

In  comparison  with  a  classical  plate  solution  (CPT),  an  element  aspect  ratio  of  10  -  which 
is  considered  in  the  regime  of  thick  plates  -  shows  an  increased  flexibility  attributed  to 
transverse  normal  and  shear  deformation  effects.  In  the  thin  plate  regime  where  a/h  >  100, 
the  maximum  deflections  are  shown  to  agree  with  the  classical  solution  for  several  reduced 
orders  of  magnitude  in  plate  thickness.  Above  a/h  =  104,  the  solution  deteriorates  due  to 
locking  which  is  clearly  evident  at  a/h  =  105.  Thus,  with  the  present  use  of  the  Pian-Tong 
hexahedral,  most  expected  layer  aspect  ratios  can  be  accurately  simulated. 

In  modeling  composite  laminates,  as  the  number  of  plies  being  modeled  in  a  single  layer 
increases,  the  effects  of  membrane-bending  and  twisting  coupling  diminish  [42].  As  these 
coupling  phenomena  vanish,  the  behavior  of  the  ply  group  tends  towards  the  behavior  of 
a  single  equivalent  homogenous  specially  orthotropic  layer.  To  demonstrate  convergence,  a 
highly  orthotropic  layer  with  properties  listed  below  is  analyzed. 


Ei 

=  40.0E6 

£2 

=  1.0E6 

Ez 

=  1.E6 

G23 

=  5.0E5 

G\z 

=  2.0E5 

G\ 2 

=  5.0E5 

V23 

=  0.25 

Viz 

=  0.25 

Vll 

=  0.25 

Plate  dimensions  were  selected  as  10  x  10  x  0.05,  and  a  uniform  normal  pressure  was  applied. 
Comparison  to  the  exact  solution  yields  a  rapidly  convergent  solution,  as  shown  in  Table  8. 

Table  8.  Deflections  of  an  equivalent  homogeneous  orthotropic  plate. 


Model 

W Max  /W Exact 

2  x  2f 

0.7911 

8  x  8f 

0.9968 

10  x  lOf 

0.9981 

16  x  16f 

1.0002 

To  assess  the  behavior  of  the  present  element  in  the  case  of  significant  coupling  between 
membrane  and  bending  deformations,  a  cross  ply  [0/90]  ply  group  is  analyzed.  The  finite 
element  model  was  composed  of  two  layers  to  assign  the  different  ply  properties  to  each 
layer.  Using  the  orthotropic  material  with  a  material  aspect  ratio  of  40  and  analyzing  a 
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thin  rectangular  plate  with  a/b  =  2  and  h  =  0.01  under  a  uniform  load,  Table  9  shows  the 
convergence  of  the  solution  compared  to  a  classical  solution  presented  by  Jones  [42], 

Table  9.  Deflections  of  a  cross-ply  [0/90]  laminate. 


Model 

Wmcix/W Exact 

2  x  2f 

0.8220 

4  x  4f 

0.9601 

8  x  8f 

0.9598 

10  x  lOf 

0.9599 

16  x  16f 

0.9604 

The  results  show  that  the  three-dimensional  formulation  converges  rapidly  to  a  maximum 
deflection  4%  less  than  the  series  solution  derived  utilizing  two-dimensional  plane  stress  as¬ 
sumptions  for  the  maximum  center  deflections  in  a  highly  coupled,  cross-ply  [0/90]  ply  group. 
For  multiple  cross-ply  or  angle-ply  groups  modeled  in  a  single  effective  layer,  the  combined 
ply  properties  exhibit  diminishing  coupling  and  behave  increasingly  as  a  homogeneous  or¬ 
thotropic  layer. 


The  essential  behavior  of  the  solid  continuum  element  incorporated  into  RESTRAN  has 
been  shown  in  a  variety  of  plate  configurations  which  have  investigated  solution  accuracy  as 
a  function  of  element  aspect  ratio,  material  orthotropy,  coupling  effects,  and  the  prediction 
of  instability  buckling.  These  results  have  been  quantified  through  comparison  with  exact 
solutions. 

The  following  illustration  in  modeling  composite  sublaminates  is  presented  without  com¬ 
parison  to  any  known  exact  solution.  A  rectangular  laminated  plate  with  a  central  el¬ 
liptical  delamination  is  selected  to  demonstrate  convergence  and  contact  constraints  in  a 
more  complicated  case  of  instability.  All  models  were  generated  using  the  RESTRAN 
♦MODEL  GENERATION  feature  with  the  parameter  NDIV  used  to  determine  the  level 
of  discretization.  The  laminate  contains  96  plies  of  nominal  thickness  0.0052  in  with  a 
near  surface  delamination  between  the  eighth  and  ninth  plies.  The  layup  is  given  by 
[±45/06||(±45)5/9012/(±45)5/06/  ±  45/(±45/06/(±45)5/906)s],  where  ||  indicates  the  lo¬ 
cation  of  the  delamination.  The  plate  dimensions  are  15  in  x  15  in  x  0.4992  in,  and  the 
elliptical  delamination  is  centrally  located  with  semimajor  axis  given  by  a  =  10.0  in  and 
semiminor  axis  by  b  =  7.5  in.  Undamaged  ply  properties  are  given  by 

£j  =  181. 0E9  E2  =  10.3E9  E3  =  10.3E9 

G23  =  7.17E9  Gn  =  1.0E9  Gn  =  1.0E9 

i/23  =  0.28  1^13  =  0.28  ^12  =  0.28 

Reduced,  isotropic  properties  were  assigned  to  the  elements  above  the  delamination  to  sim¬ 
ulate  impact  damage.  These  properties  are  given  by 
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E  =  1.0E9  G  =  3.846E8  v  -  0.3 


The  entire  plate  model  is  shown  in  Figure  35. 


Figure  35.  Laminated  plate  with  elliptical  delamination. 

The  full  model  is  reduced  to  a  half  model  by  assuming  that  the  buckling  modes  will  have  the 
x-axis  as  a  plane  of  symmetry.  The  tolerance  for  assessing  relative  motion  of  the  opposing 
delamination  surfaces  was  set  to  zero.  Figure  36  shows  the  convergence  ol  the  unconstrained 
and  constrained  buckling  modes  and  the  associated  eigenvalues. 

Buckling  Mode 


Figure  36(a).  NDIV  =  4.  No  constraints  applied.  A  =  2.0525^5. 
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Figure  36(b).  NDIV  =  4.  Compatibilty  constraints  enforced.  A  =  3.1637^5. 


Figure  36(c).  NDIV  =  8.  No  constraints  applied.  A  =  1.4667E5. 


Figure  36(d).  NDIV  =  8.  Compatibility  constraints  enforced.  A  =  2.1523E5. 


Figure  36(e).  NDIV  =  12.  No  constraints  enforced.  A  =  1.3864FI5. 


Figure  36(g).  NDIV  —  16.  No  constraints  applied.  A  —  1 . 395 6^F5 . 


Figure  36(h).  NDIV  =  16.  Compatibility  constraints  applied.  A  =  2.0527E5. 


The  eigenvalues  show  a  rapid  convergence  and  are  tabulated  in  Table  10.  In  a  linear  analysis, 
the  mode  shapes  are  of  arbitrary  sign  and  magnitude  and  can,  thus,  alternate  between  equiv¬ 
alent  configurations.  The  minimum  energy  state  corresponding  to  the  fundamental  buckling 
mode  is  a  function  of  geometry,  loading,  and  material  properties.  For  plate  geometries,  the 
buckling  mode  can  be  described  by  the  number  of  half  waves  lengths  exhibited  by  the  mode 
shape  along  each  inplane  coordinate,  (n,m).  Depending  on  the  value  of  the  independent 
plate  variables,  different  modal  buckling  patterns  can  be  elicited.  In  a  finite  element  analy¬ 
sis  the  degree  of  discretization  becomes  a  variable  in  resolving  the  mode  shape,  and  in  the 
example,  for  the  first  three  levels  of  model  refinement,  the  mode  shapes  demonstrate  a  (2,1) 
pattern.  At  the  highest  level  of  discretization  however,  a  mode  crossover  is  seen  in  Figure 
36(g),  where  the  mode  shape  changes  from  a  (2,1)  to  a  (3,1)  pattern  at  essentially  the  same 
critical  load  as  that  determined  in  the  model  shown  in  Figure  36(e). 

Table  10.  Convergence  of  unconstrained  and  constrained  critical  buckling  loads. 


NDIV 

Unconstrained 

Constrained 

4 

2.0525E+5 

3.1637E+5 

8 

1.4667E+5 

2.1523E+5 

12 

1.3864E+5 

1.9803E+5 

16 

1.3956E+5 

2.0527E+5 

The  above  set  of  standard  deformation  and  buckling  problems  demonstrate  an  accurate 
simulation  of  plate-like  bending  and  instability  behavior  using  the  incorporated  hexahedral 
element  while  calculating  three-dimensional  stress  states  for  material  failure  to  be  predicted 
using  the  same  model. 
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10  Practical  Analysis  of  Residual  Strength 

Because  of  the  complex  state  of  internal  damage  in  impacted  composites,  the  practical  ap¬ 
plication  of  any  residual  strength  predictive  methodology  requires  a  simplified  definition  of 
the  type  of  damage  and  its  spatial  distribution.  This  definition  is  required  to  analyze  ex¬ 
isting  structures  exposed  to  impact  and  new  design  concepts  with  assumed  internal  damage 
states.  The  accuracy  in  predicting  residual  strength  is  entirely  dependent  on  the  accuracy  in 
which  internal  damage  can  be  resolved.  For  in-service  structural  components,  a  depot-level 
inspection  must  be  capable  of  discerning  interal  damage  accurately  enough  to  safely  assess 
repair/no-repair  options  for  continued  service  applications.  For  design,  specific  characteri¬ 
zations  of  maximum  sustained  impact  damage  need  to  be  developed  to  assess  survivability 
and  damage  tolerance  of  structures  to  expected  impact  threats. 


10.1  Damage  Characterization 

For  a  given  impact  event,  the  extent  of  impact  damage  can  be  determined  from  a  simulation  of 
the  impact  event  based  on  first  principles  in  elastodynamics,  nonlinear  material  constitutive 
relations,  and  contact  dynamics.  A  basic  measure  is  the  amount  of  energy  absorbed  due  to 
impact.  This  may  be  expressed  as 

u = v- 1!  fm + »  jC  Ftdt  -  h  (I! Fdt ) 2  (l63) 

where  v0  is  the  velocity  of  the  impactor,  M  is  the  mass  of  the  impactor,  F  is  the  impact 
force,  g  is  the  acceleration  due  to  gravity,  and  t  is  the  time  duration  of  the  impact  event.  In 
turn,  this  energy  can  be  differentiated  into  components,  Ue  and  Up >,  where  Ue  is  the  energy 
converted  into  elastic  deformation  of  the  specimen,  and  Ud  is  the  component  of  energy 
absorbed  in  irreversible  daqmage  causing  mechanisms.  Assuming  basic  damage  modes  in 
composite  materials,  the  damage  energy  may  be  further  decomposed  as 

UD  =  UM  +  UF  +  UC  (164) 

where  Um  is  the  energy  absorbed  in  matrix  cracking,  Up  is  the  energy  consumed  in  figer 
breakage,  and  Uc  is  the  energy  required  to  cause  the  creation  of  fracture  surfaces  in  the 
form  of  interlaminar  delaminations  and  transverse  cracks.  The  apportioning  of  energy  to 
each  damage  mode  must  be  assumed,  measured,  or  calculated  to  analytically  predict  the 
three-dimensional  spatial  distribution  of  damage  [63]. 

The  increasing  development  and  application  of  health  monitoring  technologies  such  as  em¬ 
bedded  fiber  optic  sensors  will  permit  an  increasingly  accurate  characterization  of  internal 
damage  states  due  to  sustained  impacts.  Continued  improvements  in  established  nonde¬ 
structive  evaluation  methods  such  as  acoustics/ultrasonics  [83-89],  thermography  [90,91], 
and  computed  tomography  [92-95]  are  providing  increasingly  accurate  resolutions  of  internal 
damage  states  and  the  ability  to  distinguish  different  damage  modes. 
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The  increasing  ability  to  characterize  internal  damage  states  in  existing  composite  struc¬ 
tures  exposed  to  impact  will  better  define  damage  tolerant  design  objectives  and  provide 
more  accurate  initial  damage  conditions  for  input  to  analytical  methodologies  such  as  RE- 
STRAN  for  predicting  residual  strength. 


10.2  Demonstration  Problem 

To  illustrate  the  execution  of  the  RESTRAN  analysis  program,  a  model  of  an  elastically 
supported  laminated  composite  face  plate  is  analyzed.  Due  to  the  lack  of  data  of  sufficient 
resolution  regarding  the  spatial  distribution  of  specific  internal  damage  modes  and  subse¬ 
quent  experimental  determination  of  residual  strength,  a  specific  state  of  internal  damage  is 
assumed  for  the  following  example.  The  geometry  and  applied  loading  is  shown  in  Figure 
37. 


Figure  37.  Geometry  and  loading  of  an  elastically  supported  composite  plate. 


Material  properties  were  selected  as  S2-Glass/3501  Epoxy  tape  with  a  nominal  ply  thickness 
of  0.0052  in.  The  material  elastic  properties  are  given  by 

Ex  =  7.150E6  Psi  E2  =  2.13E6  Psi  Ez 

G2Z  =  0.71E6  Psi  Gw  =  0.98E6  Psi  Gn 

v2Z  =  0.499  uxz  =  0.306  u12 

with  strengths  given  by 

XTen  =  2.43E5  Psi  XComP  =  1.77E5  Psi  YTen  =  7.0E3  Psi 

Ycomp  =  3.06E4  Psi  ZTen  =  7.03E3  Psi  ZComV  =  3.5E4  Psi 

R  =  1.7E4  Psi  S  =  1.57E4  Psi  T  =  1.57E4  Psi 


=  2.13E6  Psi 
=  0.98E6  Psi 
=  0.296 
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The  face  plate  is  composed  of  96  plies  and  is  assumed  to  have  been  subjected  to  a  normal 
impact  leading  to  an  inner  circular  region  where  two  delaminations  exist  located  eight  plies 
or  0.0416  in  from  the  outer  and  inner  laminate  surface.  Furthermore,  the  material  proper¬ 
ties  of  these  outer  two  layers  are  assumed  to  have  experienced  fiber  damage  such  that  the 
longitudinal  modulus  has  been  reduced  to  the  matrix-dominated  properties  of  the  transverse 
moduli.  The  layup  and  location  of  these  delaminations,  indicated  by  ’||’  in  the  layup  descrip¬ 
tion,  is  given  by  [±45/06  ||  ±455/906 ± 455/06 ±452/0s/ ±452/06±455/90e/ ±45s  ||  0^/ ±45]. 
A  residual  strength  analysis  is  performed,  which  is  selected  to  evaluate  both  delamination 
instability  and  material  failure.  For  this  example,  the  maximum  stress  failure  criteria  is  used 
with  a  nonspecific  damage  law  that  assigns  zero  modulus  to  plies  experiencing  failure  with¬ 
out  regard  to  failure  mode.  Abbreviated  input  and  output  files  are  shown  in  the  following 
subsections,  followed  by  a  graphical  presentation  of  the  progression  of  failure. 


10.2.1  Input  Data  File  for  a  Residual  Strength  Problem 


♦HEADING 

MODEL  OF  AN  ELASTICALLY  SUPPORTED  COMPOSITE  PLATE 
❖  ❖ 

**  [45/-45/0_6 | | (45/-45)_5/90_6/(45/-45)_5/0_6/(45/-45)_2/0_8/ 
**  (-45/45)_2/0_6/(-45/45)_5/90_6/(-45/45)_5| |0_6/-45/45] 

** 

**  LAYUP  TAPE  THICKNESS  =  0 . 0052 ’ ’ 

** 

**  TWO  INTERIOR  DELAMINATIONS  MODELLED  AT 
**  LAYER  INTERFACES  2 | 3  AND  3  1 4 
** 

❖SOLUTION,  METH  =  CMB 
6,15,1.05 
❖❖ECHO 
❖❖PREPASS 
❖❖NODE  PRINT 
**  D,M 

❖PLY  FAILURE  PRINT 
❖parameter  directory  =  current 
❖parameter  lapfail 
❖parameter  status 
❖MEMORY  ALLOCATION 
BAND 

❖GRAPHICS,  format  =  mathematica 
2,  2,  1.0 
❖  ❖ 

**  NODE  DEFINITIONS 
❖  * 

❖NODE 
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**  TOP 

FRONT 

COMPONENT 

1 

5.0 

-10.0 

10.0 

2 

5.8 

-10.0 

10.0 

3 

6.3 

-10.0 

10.0 

3511 

9.3375 

10.0 

-5.9189 

3512 

9.4375 

10.0 

-5.9189 

3513 

9.5 

10.0 

-5.9189 

** 

**  BOTTOM,  FRONT 

COMPONENT 

** 

3601 

5.0 

-10.0 

-10.0 

3602 

5.8 

-10.0 

-10.0 

3603 

6.3 

-10.0 

-10.0 

4711 

9.3375 

-8.6396 

-5.9189 

4712 

9.4375 

-8.6396 

-5.9189 

4713 

9.5 

-8.6396 

-5.9189 

♦♦ 

**  NODES,  ELEMENTS,  AND  NODE  SETS  FOR  DELAMINATED  VERTICAL 
**  COMPONENT  OBTAINED  FROM  THE  FOLLOWING  *MODEL  GENERATION 
♦  *  STATEMENT: 


** 

♦MODEL  GENERATION 

** 

8,4,10.0,10.0,4,4,4,4 

** 

2. 0,1, 1,1, 1,1 

** 

0.0,  0.0,  1.0 

** 

0.0,  1.0,  0.0 

** 

1.0,  0.0,  0.0 

** 

1.0,10.0 

** 

1,5,250,5000,5000 

** 

0.95,9.875 

** 

2,6,255,7000,7000 

** 

0.85,9.625 

** 

3,7,265,11000,11000 

** 

0.8,9.50 

** 

4,8,270,13000,13000 

♦NODE 


5001 

0 . 1000000E+02 

O.OOOOOOOE+OO 

O.OOOOOOOE+OO 

5002 

0 . 1000000E+02 

0.3750000E+00 

O.OOOOOOOE+OO 

5003 

0 . 1000000E+02 

0.3750000E+00 

0.3750000E+00 
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13519  0.9500000E+01  0. 1000000E+02 

13520  0.9500000E+01  0. 1000000E+02 

13521  0 . 9500000E+01  0. 1000000E+02 

** 

**  DEFINE  GENERATOR  ELEMENTS 


0.7279306E+01 
0.8639653E+01 
0 . 1000000E+02 


** 

* ELEMENT,  LAYUP  =  7,  elset  =  300 
** 

**  top  front  elastic  support 
** 


1, 

401, 

402, 

902,  901, 

1, 

2, 

102, 

101 

101, 

501, 

502, 

802,  801, 

401, 

402, 

902, 

901 

201, 

601, 

602, 

702,  701, 

501, 

502, 

802, 

801 

301, 

901, 

902, 

1002,  1001, 

101, 

102, 

202, 

201 

401, 

** 

**  top 

j|ca|c 

1001, 

1002, 

1102,  1101, 

201, 

202, 

302, 

301 

back 

elastic 

support 

501, 

1601, 

1602, 

1702,  1701, 

1201, 

1202, 

1302, 

1301 

601, 

1701, 

1702, 

1802,  1801, 

1301, 

1302, 

1402, 

1401 

701, 

1801, 

1802, 

1902,  1901, 

1401, 

1402, 

1502, 

1501 

801, 

2001, 

2002, 

2102,  2101, 

1801, 

1802, 

1902, 

1901 

901, 

2201, 

2202, 

2302,  2301, 

2001, 

2002, 

2102, 

2101 

** 

**  bottom  front  elastic  support 


** 

1001, 

3601, 

3602, 

3702, 

3701, 

4001, 

4002, 

4102, 

4101 

1101, 

4001, 

4002, 

4102, 

4101, 

4401, 

4402, 

4502, 

4501 

1201, 

4401, 

4402, 

4502, 

4501, 

4601, 

4602, 

4702, 

4701 

1301, 

3701, 

3702, 

3802, 

3801, 

4101, 

4102, 

4202, 

4201 

1401, 

3801, 

3802, 

3902, 

3901, 

4201, 

4202, 

4302, 

4301 

** 

**  bottom  back  elastic  support 
** 


1501, 

2401, 

2402, 

2502, 

2501, 

2801, 

2802, 

2902, 

2901 

1601, 

2501, 

2502, 

2602, 

2601, 

2901, 

2902, 

3002, 

3001 

1701, 

2601, 

2602, 

2702, 

2701, 

3001, 

3002, 

3102, 

3101 

1801, 

3001, 

3002, 

3102, 

3101, 

3201, 

3202, 

3302, 

3301 

1901, 

3201, 

3202, 

3302, 

3301, 

3401, 

3402, 

3502, 

3501 

** 


**  ELEMENT  GENERATION  FOR  ELASTIC  SUPPORTS 
** 


*ELGEN 

1, 

101, 


,  ELID  =  300 


12,1,1 

12,1,1 


j  >  y  9  9  y 

>  >  y  9  9  9 
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201, 


12,1,1 


9  9  9  9  9  9 


1701,  12,1,1,,,,,, 

1801,  12,1,1,,,,,, 

1901,  12,1,1,,,,,, 

** 

**  CREATE  NODE  SETS  FOR  EQUIVALENCING 
♦♦ 


♦NSET, 

NSID  = 

10 

13, 

113, 

213, 

313, 

413, 

913, 

1013, 

1113 

513, 

813, 

613, 

713 

♦NSET, 

NSID  = 

20 

13300, 

13319, 

13338, 

13357, 

13299, 

13318, 

13337, 

13356 

13298, 

13317, 

13297, 

13316 

♦NSET, 

NSID  = 

30 

1213, 

1313, 

1413, 

1513, 

1613, 

1713, 

1813, 

1913 

2013, 

2113, 

2213, 

2313 

♦NSET, 

NSID  = 

40 

13464, 

13483, 

13502, 

13521, 

13463, 

13482, 

13501, 

13520 

13500, 

13519, 

13499, 

13518 

♦NSET, 

NSID  = 

50 

4613, 

4713, 

4413, 

4513, 

4013, 

4113, 

4213, 

4313 

3613, 

3713, 

3813, 

3913 

♦NSET, 

NSID  = 

60 

13285, 

13304, 

13284, 

13303, 

13283, 

13302, 

13321, 

13340 

13282, 

13301, 

13320, 

13339 

♦NSET, 

NSID  = 

70 

3413, 

3513, 

3213, 

3313, 

2813, 

2913, 

3013, 

3113 

2413, 

2513, 

2613, 

2713 

♦NSET, 

NSID  = 

80 

13487, 

13506, 

13486, 

13505, 

13447, 

13466, 

13485, 

13504 

13446, 

13465, 

13484, 

13503 

**  EQUIVALENCE  NODE  SETS 
♦♦ 

♦EQUIVALENCE 
10,  20 
30,  40 
50,  60 
70,  80 
♦♦ 

**  ELEMENTS  IN  VERTICAL  COMPONENT  FROM  PRIOR  ♦MODEL  GENERATION  RUN 
♦♦ 

♦ELEMENT,  LAYUP  =  1  ORIENTATION  =  1 

5001  7001  7004  7003  7002  5001  5004  5003 


5002 
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5002 

7006 

7005 

7004 

7001 

5006 

5005 

5004 

5001 

5003 

7007 

7006 

7001 

7008 

5007 

5006 

5001 

5008 

5126 

7111 

7112 

7144 

7143 

5111 

5112 

5144 

5143 

5127 

7112 

7113 

7145 

7144 

5112 

5113 

5145 

5144 

5128 

7113 

7082 

7114 

7145 

5113 

5082 

5114 

5145 

♦ELEMENT, 

LAYUP  = 

4  ORIENTATION  = 

1 

5129 

7114 

7115 

7147 

7146 

5114 

5115 

5147 

5146 

5130 

7115 

7116 

7148 

7147 

5115 

5116 

5148 

5147 

5131 

7116 

7117 

7149 

7148 

5116 

5117 

5149 

5148 

5482 

7499 

7500 

7519 

7518 

5499 

5500 

5519 

5518 

5483 

7500 

7501 

7520 

7519 

5500 

5501 

5520 

5519 

5484 

7501 

7502 

7521 

7520 

5501 

5502 

5521 

5520 

♦ELEMENT, 

LAYUP  = 

2  ORIENTATION 

= 

1 

7001 

11001 

11004 

11003 

11002 

7001 

7004 

7003 

7002 

7002 

11006 

11005 

11004 

11001 

7006 

7005 

7004 

7001 

7003 

11007 

11006 

11001 

11008 

7007 

7006 

7001 

7008 

7126 

11111 

11112 

11144 

11143 

7111 

7112 

7144 

7143 

7127 

11112 

11113 

11145 

11144 

7112 

7113 

7145 

7144 

7128 

11113 

11082 

11114 

11145 

7113 

7082 

7114 

7145 

♦ELEMENT, 

LAYUP  = 

5  ORIENTATION  = 

1 

7129 

11114 

11115 

11147 

11146 

7114 

7115 

7147 

7146 

7130 

11115 

11116 

11148 

11147 

7115 

7116 

7148 

7147 

7131 

11116 

11117 

11149 

11148 

7116 

7117 

7149 

7148 

7482 

11499 

11500 

11519 

11518 

7499 

7500 

7519 

7518 

7483 

11500 

11501 

11520 

11519 

7500 

7501 

7520 

7519 

7484 

11501 

11502 

11521 

11520 

7501 

7502 

7521 

7520 

♦ELEMENT , 

LAYUP  = 

3  ORIENTATION 

= 

1 

11001 

13001 

13004 

13003 

13002 

11001 

11004 

11003 

11002 

11002 

13006 

13005 

13004 

13001 

11006 

11005 

11004 

11001 

11003 

13007 

13006 

13001 

13008 

11007 

11006 

11001 

11008 

11126 

13111 

13112 

13144 

13143 

11111 

11112 

11144 

11143 

11127 

13112 

13113 

13145 

13144 

11112 

11113 

11145 

11144 
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11128 

13113 

13082 

13114 

13145 

11113 

11082 

11114 

11145 

♦ELEMENT, 

LAYUP  = 

6  ORIENTATION 

= 

1 

11129 

13114 

13115 

13147 

13146 

11114 

11115 

11147 

11146 

11130 

13115 

13116 

13148 

13147 

11115 

11116 

11148 

11147 

11131 

13116 

13117 

13149 

13148 

11116 

11117 

11149 

11148 

11482 

13499 

13500 

13519 

13518 

11499 

11500 

11519 

11518 

11483 

13500 

13501 

13520 

13519 

11500 

11501 

11520 

11519 

11484 

13501 

13502 

13521 

13520 

11501 

11502 

11521 

11520 

*♦ 

**  RELAX  TOLERANCE  OF  DEFORMED  GEOMETRY  CHECKS 
** 

♦DEFORMED  GEOMETRY 

20.0 


♦♦ 


**  LAMINATE 

** 

DESCRIPTION 

♦LAYER,  LAYUP  =  1 

3,  0.0052, 

45.0 

3,  0.0052, 

-45.0 

3,  0.0312, 

0.0 

♦LAYER,  LAYUP  =  2 

1,  0.0260, 

45.0 

1,  0.0260, 

-45.0 

1,  0.0312, 

90.0 

1,  0.0260, 

45.0 

1,  0.0260, 

-45.0 

1,  0.0312, 

0.0 

1,  0.0104, 

-45.0 

1,  0.0104, 

45.0 

1,  0.0416, 

0.0 

1,  0.0104, 

-45.0 

1,  0.0104, 

45.0 

1,  0.0312, 

0.0 

1,  0.0260, 

-45.0 

1,  0.0260, 

45.0 

1,  0.0312, 

90.0 

1,  0.0260, 

-45.0 

1,  0.0260, 

45.0 

♦LAYER,  LAYUP=3 

3,  0.0312, 

0.0 

3,  0.0052, 

-45.0 

3,  0.0052, 

45.0 

♦LAYER,  LAYUP  =  4 

1,  0.0052, 

45.0 
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1,  0.0052, 

-45.0 

1,  0.0312, 

0.0 

♦LAYER,  LAYUP  =  5 

1,  0.0260, 

45.0 

1,  0.0260, 

-45.0 

1,  0.0312, 

90.0 

1,  0.0260, 

45.0 

1,  0.0260, 

-45.0 

1,  0.0312, 

0.0 

1,  0.0104, 

-45.0 

1,  0.0104, 

45.0 

1,  0.0416, 

0.0 

1,  0.0104, 

-45.0 

1,  0.0104, 

45.0 

1,  0.0312, 

0.0 

1,  0.0260, 

-45.0 

1,  0.0260, 

45.0 

1,  0.0312, 

90.0 

1,  0.0260, 

-45.0 

1,  0.0260, 

45.0 

♦LAYER,  LAYUP=6 

1,  0.0312, 

0.0 

1,  0.0052, 

-45.0 

1,  0.0052, 

45.0 

♦LAYER,  LAYUP  =  7 
2,  1.0,  0.0 
♦  ♦ 

♦♦  MATERIAL  DEFINITIONS 
♦  ♦ 

**  COMPOSITE  PLY  PROPERTIES  (S2-GLASS/3501  EPOXY) 
** 

♦MATERIAL,  MATID  =  1 

7.150E6,  2.13E6,  2.13E6,  0.98E6,  0.71E6,  0.98E6 
0.306,  0.499,  0.296 
♦FAILURE  CRITERIA,  FCID  =  1 
MAX-STRESS 

2.43E5,  1.77E5,  7.0E3,  3.06E4,  7.0E3,  3.5E4 
1.7E4,  1.57E4,  1.57E4 
♦DAMAGE  LAW,  DLID  =  1 
NULL 
♦  ♦ 

♦♦  REDUCED  PROPERTIES  WITHIN  INNER  REGION 
♦  ♦ 

♦MATERIAL,  MATID  =  3 

2.13E6,  2.13E6,  2.13E6,  8.15E5,  8.15E5,  8.15E5 
0.306,  0.306,  0.306 
♦FAILURE  CRITERIA,  FCID  =  3 
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MAX-STRESS 

2.43E5,  1.77E5,  7.0E3,  3.06E4,  7.0E3,  3.5E4 
1.7E4,  1.57E4,  1.57E4 
♦DAMAGE  LAW,  DLID  =  3 
NULL 
** 

**  METALLIC  PROPERTIES 
** 

♦MATERIAL,  MATID  =  2 

1.0E8,  1.0E8,  1.0E8,  0.3846E8,  0.3846E8,  0.3846E8 
0.3,  0.3,  0.3 

♦FAILURE  CRITERIA,  FCID  =  2 
MAX-STRESS 

5.3E5,  2.2E6,  5.3E5,  2.2E6,  5.3E5,  2.2E6 
4.8E5,  4.8E5,  4.8E5 
♦DAMAGE  LAW,  DLID  =  2 
NULL 
** 

♦♦  ESTABLISH  LOCAL  COORDINATE  SYSTEM 
** 

♦ORIENTATION 

1,  0.0,  1.0,  0.0,  0.0,  0.0,  1.0 

** 

♦♦  NODE  SET  DEFINITION  FOR  BOUNDARY  CONSTRAINT  INPUT 
** 

♦NSET,  NS ID  =100 

1,  101,  201,  301,  401,  501,  601,  701,  801 

901,  1001,  1101 
♦NSET,  NSID  =110 

1201,  1301,  1401,  1501,  1601,  1701,  1801,  1901,  2001 
2101,  2201,  2301 
♦NSET,  NSID  =120 

2401,  2501,  2601,  2701,  2801,  2901,  3001,  3101,  3201 
3301,  3401,  3501 
♦NSET,  NSID  =130 

3601,  3701,  3801,  3901,  4001,  4101,  4201,  4301,  4401 
4501,  4601,  4701 
♦BOUND ARY2 
100,  1,  1 
110,  1,  2 
120,  1,  3 
130,  1,  3 

** 

♦♦  INPUT  NODE  SETS  USED  TO  DEFINE  DELAMINATION  PLANES 
** 

♦NSET,  NSID  =  955 

7001  7002  7003  . . .  7143  7144  7145 
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*NSET,  NS ID  =  960 

11001  11002  11003  ...  11143  11144  11145 

** 

**  INPUT  DELAMINATIONS  AS  SPECIFIED  NODE  SETS 
** 

*DELAMINATION 
955  960 
** 

**  EXCLUDE  ELEMENTS  COMPRISING  ELASTIC  SUPPORTS 
**  FROM  FAILURE  PREDICTION 
** 

♦EXCLUDE  ELEMENT 
300 
** 

**  biaxial  y-z  plane  loading 
** 

♦CLOAD 

** 

**  normal  z-direction  loads 
** 

5300,  3,  -. 1000000E+02 
5319,  3,  -.1000000E+02 
5338,  3,  -. 1000000E+02 


13465,  3,  0 . 1000000E+02 

13484,  3,  0 . 1000000E+02 

13503,  3,  0 . 1000000E+02 

** 

**  transverse  y-direction  loads 
** 

5282,  2,  0.2000000E+02 

5283,  2,  0.2000000E+02 

5284,  2,  0.2000000E+02 


13519,  2,  -.2000000E+02 

13520,  2,  -.2000000E+02 

13521,  2,  - . 2000000E+02 

♦ENDDATA 


10.2.2  Output  Data  File  for  a  Residual  Strength  Problem 


®®®@@®®@@®@@@®@®@®@@@@@®@@@®®@®@®@®@®@@@®@@@®@ 


m  @@ 

@®  U.S.  ARMY  RESEARCH  LABORATORY  ®® 

®@  m 

@0  RESTRAN  ®® 

®®  @® 

@@  RESIDUAL  STRENGTH  ANALYSIS  OF  IMPACT  ®@ 
®®  DAMAGED  COMPOSITE  LAMINATES  ®@ 

@®  ®® 

®®  VERSION  1.0  ®® 

®®  ®® 


<§@@@@@@@<a@@@@@@@@®®@®@®®@@@@@®®@®@@@@®@@@®®®@@ 

***  MESSAGE:  ELEMENT  NODE  ORDER  IS  BEING  CONVERTED  TO  RESTRAN  FORMAT 

***  WARNING:  MATERIAL  ID  1  HAS  A  MIXED  MODE  FAILURE  CRITERIA 

ASSOCIATED  WITH  A  SINGLE  MODE  DAMAGE  LAW  (NULL) . 

NULL  ACCEPTED 

***  WARNING:  MATERIAL  ID  3  HAS  A  MIXED  MODE  FAILURE  CRITERIA 

ASSOCIATED  WITH  A  SINGLE  MODE  DAMAGE  LAW  (NULL) . 

NULL  ACCEPTED 

***  WARNING:  MATERIAL  ID  2  HAS  A  MIXED  MODE  FAILURE  CRITERIA 

ASSOCIATED  WITH  A  SINGLE  MODE  DAMAGE  LAW  (NULL)  . 

NULL  ACCEPTED 


###################################### 
###  ### 

###  MODEL  DEFINITION### 
###  ### 

###################################### 


FORCE/MOMENT  RESULTANTS  AT  ORIGIN: 

I R |  =  0 . 400E+02  IMI  =  0.456E+03 
Rx  =  0 . OOOE+OO  Ry  =  O.OOOE+OO  Rz.=  -.400E+02 
Mx  =  0 . 237E+03  My  =  -.390E+03  Mz  =  O.OOOE+OO 
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Xo  =  0.975E+01  Yo  =  0.592E+01  Zo  =  O.OOOE+OO 


***  END  OF  MODEL  DEFINITION  *** 


MODEL  SIZE  PARAMETERS 


NUMBER  OF  ELEMENTS  =  1692 

NUMBER  OF  NODES  =  2950 

DEGREES  OF  FREEDOM  =  8850 

SYSTEM  BANDWIDTH  =  2892 

############################################### 
##  ## 

##  BEGIN  FAILURE  ANALYSIS## 
##  ## 


############################################### 


@@  @@ 
@@  PRIMARY  ANALYSIS  CYCLE.  PASS  NUMBER  1  @<§ 


«<  ANALYZING  DELAMINATION  SET:  1  »> 

.«<  CONTACT  ITERATION  NUMBER:  1  »> 

CONVERGED  EIGENVALUES: 

NUMBER  LAMBDA 

1  0 . 1825575E+02 

2  0 . 1482253E+02 

3  0 . 1040733E+02 

4  0 . 6461659E+01 

DELAMINATION  1  AT  ITERATION  1  IS  EXHIBITING 
A  CONTACT  COMPATIBILITY  OF  100.00  PERCENT 

«<  ANALYZING  DELAMINATION  SET:  2  »> 

«<  CONTACT  ITERATION  NUMBER:  1  »> 

CONVERGED  EIGENVALUES: 

NUMBER  LAMBDA 
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1  0 . 1826281E+02 

2  0 . 1488730E+02 

3  0. 1037982E+02 

4  0.6463456E+01 

DELAMINATION  2  AT  ITERATION  1  IS  EXHIBITING 
A  CONTACT  COMPATIBILITY  OF  100.00  PERCENT 


* ********************************** ************** 


**  ** 

**  ALGORITHMIC  PATH  FOR  MINIMUM  LOAD  INCREMENT  ** 

**  TO  NEXT  FAILURE:  DELAMINATION  INSTABILITY  ** 

**  ** 

**  MATERIAL  BUCKLING  ** 

**  -  -  ** 

**  SCALE:  0.540E+02  0.646E+01  ** 

**  ** 


j****************:****:*:***#:*****:****:*************** 

******************************************* 

**  ** 

**  SUBLAMINATE  BUCKLING  FAILURE  ANALYSIS  ** 

**  ** 

******************************************* 


*  DELAMINATION  DEFINED  BY  NODE  SET  955 
IS  PREDICTED  TO  BUCKLE  AT  AN  APPLIED 
LOAD  FACTOR  OF  0.646E+01 

*  FAILURE  IS  PREDICTED  TO  INCLUDE  ALL  ELEMENTS  ALONG 
POSITIVE  NORMAL  TO  THE  DELAMINATION  PLANE. 

*  THE  FOLLOWING  ELEMENTS  HAVE  BEEN  DEGRADED 
OR  FAILED  VIA  LOCAL  SUBLAMINATE  BUCKLING: 

5001  5002  5003  5004  . . .  5157  5158  5159 

*  LOCAL  SUBLAMINATE  BUCKLING  IS  PREDICTED. 

ELEMENT  PLY  FAILURE  STATUS  AT  CYCLE  NO  1 


ELEMENT 

PLY  FAILURE  MODES: 

ID 

‘/.FIBER 

•/.MATRIX 

•/.BUCKLING 

•/.TOTAL 

5001 

0 

0 

100 

100 

5002 

0 

0 

100 

100 

5003 

0 

0 

100 

100 

5160 
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5158 

0 

0 

49 

49 

5159 

0 

0 

49 

49 

5160 

0 

0 

49 

49 

*  THE  FOLLOWING  DELAMINATION  NODE  SET(S)  WILL 
BE  REMOVED  DUE  TO  BUCKLING  FAILURE: 

955 


339  DEGREES  OF  FREEDOM  ARE  CURRENTLY  ASSOCIATED  WITH 
ZERO  STIFFNESS  AND  HAVE  BEEN  REMOVED 

«<  ANALYZING  DELAMINATION  SET:  1  »> 

«<  CONTACT  ITERATION  NUMBER:  1  »> 

CONVERGED  EIGENVALUES: 

NUMBER  LAMBDA 

1  0.2931583E+01 

2  0 . 2453289E+01 

3  0 . 1653048E+01 

4  0 . 1043609E+01 

DELAMINATION  1  AT  ITERATION  1  IS  EXHIBITING 
A  CONTACT  COMPATIBILITY  OF  100.00  PERCENT 


*  *  *  *  *  *  ** * *  *  **  *  **  *  *  *  *********************  *  *  **  *  *****  * 


**  ** 

**  ALGORITHMIC  PATH  FOR  MINIMUM  LOAD  INCREMENT  ** 

**  TO  NEXT  FAILURE:  DELAMINATION  INSTABILITY  ** 

**  ** 

**  MATERIAL  BUCKLING  ** 

**  -  -  ** 

**  SCALE:  0 . 536E+02  0.674E+01  ** 

**  ** 


************************************************* 
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**  ** 

**  SUBLAMINATE  BUCKLING  FAILURE  ANALYSIS  ** 

**  ** 

*  DELAMINATION  DEFINED  BY  NODE  SET  960 
IS  PREDICTED  TO  BUCKLE  AT  AN  APPLIED 
LOAD  FACTOR  OF  0.674E+01 

*  FAILURE  IS  PREDICTED  TO  INCLUDE  ALL  ELEMENTS  ALONG 
NEGATIVE  NORMAL  TO  THE  DELAMINATION  PLANE. 

*  THE  FOLLOWING  ELEMENTS  HAVE  BEEN  DEGRADED 
OR  FAILED  VIA  LOCAL  SUBLAMINATE  BUCKLING: 

11001  11002  11003  11004  ...  11157  11158  11159  11160 

*  LOCAL  SUBLAMINATE  BUCKLING  IS  PREDICTED. 


ELEMENT  PLY  FAILURE  STATUS  AT  CYCLE  NO  2 


ELEMENT 

PLY  FAILURE  MODES: 

ID 

’/.FIBER 

'/.MATRIX 

'/.BUCKLING 

'/.TOTAL 

5001 

0 

0 

100 

100 

5002 

0 

0 

100 

100 

5003 

0 

0 

100 

100 

11158 

0 

0 

49 

49 

11159 

0 

0 

49 

49 

11160 

0 

0 

49 

49 

*  THE  FOLLOWING  DELAMINATION  NODE  SET(S)  WILL 
BE  REMOVED  DUE  TO  BUCKLING  FAILURE: 


960 


3  @0 


PRIMARY  ANALYSIS  CYCLE.  PASS  NUMBER 


«<  ANALYZING  DELAMINATION  SET:  1  »> 

«<  CONTACT  ITERATION  NUMBER:  1  »> 

CONVERGED  EIGENVALUES: 

NUMBER  LAMBDA 

1  0 . 2200270E+02 

2  0 . 1958573E+02 

3  0 . 1437475E+02 

4  0 . 8520247E+01 


************************************************* 


** 


** 


**  ALGORITHMIC  PATH  FOR  MINIMUM  LOAD  INCREMENT  ** 
**  TO  NEXT  FAILURE:  MATERIAL  DEGRADATION  ** 


** 

** 

** 

**  SCALE: 
** 


MATERIAL 
0 . 535E+02 


BUCKLING 
0 . 575E+02 


** 

** 

** 

** 

** 


******* ******************************* *********** 


********************************* 


*  * 

*  MATERIAL  FAILURE  ANALYSIS  * 

*  * 

*  ITERATION  NO.  SCALE  FACTOR  * 

*  -  - * 

*  1  0.535E+02  * 

*  * 


********************************* 


*  NUMBER  OF  ELEMENTS  DEGRADED  =  0 


*** * **************************** * 


* 

* 

* 

* 

* 

* 

* 


* 

MATERIAL  FAILURE  ANALYSIS  * 

* 

ITERATION  NO.  SCALE  FACTOR  * 
-  - * 

2  0 . 535E+02  * 

* 


********************************* 


*  NUMBER  OF  ELEMENTS  DEGRADED  =  8 
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**** ********************* ******** 


*  * 

*  MATERIAL  FAILURE  ANALYSIS  * 

*  * 

*  ITERATION  NO.  SCALE  FACTOR  * 

*  -  - * 

*  3  0 . 535E+02  * 

*  * 


* ** ** ** * **** * * * ** * ** * ** * ******** * 

*  NUMBER  OF  ELEMENTS  DEGRADED  =  30 

***  **  *  **  ****  *  *******************  * 


*  * 

♦  MATERIAL  FAILURE  ANALYSIS  * 

*  * 

*  ITERATION  NO.  SCALE  FACTOR  * 

*  -  - * 

*  4  0.535E+02  * 

*  * 


**** ******* ********* ************* 

*  NUMBER  OF  ELEMENTS  DEGRADED  =  111 

********************************* 


*  * 

*  MATERIAL  FAILURE  ANALYSIS  * 

*  * 

*  ITERATION  NO.  SCALE  FACTOR  * 

*  -  - * 

*  5  0 . 535E+02  * 

*  * 


********************************* 

*  NUMBER  OF  ELEMENTS  DEGRADED  =  211 

*  ELEMENT  FAILURE  HAS  ALTERED  MODEL  STABILITY  SUCH 
THAT  RIGID  BODY  MODES  HAVE  BEEN  DETECTED. 

TOTAL  FAILURE  IS  ASSUMED. 

ELEMENT  PLY  FAILURE  STATUS  AT  CYCLE  NO  3 


ELEMENT  PLY  FAILURE  MODES: 

ID  ‘/.FIBER  ‘/.MATRIX  ‘/.BUCKLING  ‘/.TOTAL 
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5001 

5002 

5003 


0  0  100  100 
0  0  100  100 
0  0  100  100 


11462 

11463 

11464 


0  0 
0  0 
0  0 


0  33 

0  100 

0  100 


############################################### 
##  ## 

##  ANALYSIS  PREDICTS  CATASTROPHIC  FAILURE  OF  ## 

##  THE  MODEL  AT  AN  ULTIMATE  LOAD  GIVEN  BY:  ## 

##  ## 

##  P(ULT)  =  (0.53468E+02)  *  P (INITIAL)  ## 

##  ## 

############################################### 


10.2.3  Graphical  Output  Using  MATHEMATICA 

Preface  graphic  pages  are  generated  to  show  the  color  codes  used  to  display  various  element 
failure  modes.  In  the  specific  demonstration  problem  analyzed,  complete  element  failure  is 
indicated  by  a  wireframe  depiction,  and  because  no  specific  failure  modes  were  predicted 
using  the  maximum  stress  criterion,  material  failure  is  shown  as  using  the  general  material 
failure  color  coding  scheme. 
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Both  delaminations  are  analyzed  separately  for  critical  buckling  load.  These  are  compared 
to  the  load  required  to  cause  first  ply  failure,  and  it  is  determined  that  the  outer  surface 
delamination  is  the  first  failure  to  occur.  Material  properties  for  the  elements  involved  in 
the  sublaminate  buckling  are  set  to  zero  such  that  they  appear  in  a  wireframe  depiction  in 
the  failure  state  graphic.  This  is  followed  by  a  magnified  view  of  the  buckling  mode  shown 
in  the  following  illustration. 
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ANALYSIS  CYCLE  1 

Buckling  Failure  at  Scale  =  0.646E+01 


Failure  State 


Suckling  Mode 


v 


The  next  predicted  failure  is  buckling  of  the  layer  on  the  inner  plate  surface,  as  shown 
next. 
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The  final  sequence  of  failure  events  as  fully  described  in  the  output  file  involve  a  cascade 
of  material  ply  failures.  These  occured  during  the  third  global  analysis  cycle  which,  after 
five  iterations,  had  reduced  the  overall  stiffness  such  that  rigid  body  modes  were  detected 
and  the  analysis  terminated.  The  final  failure  state  is  shown  next,  in  which  total  element 
failures  essentially  created  an  opening  through  the  inner  region  of  the  face  plate.  Mediating 
element  failures  are  shown  to  have  progressed  to  the  outer  plate  boundary.  Failure  in  the 
line  of  elements  comprising  the  outer  boundary  have  been  precluded  using  the  LAPFAIL 
parameter  which  prevents  elements  to  which  applied  external  loads  have  been  applied  from 
exhibiting  failure. 


ANALYSIS  CYCLE  3 

Material  Failure  at  Scale  =  0.535E+02 


The  final  residual  strength  prediction  yields  the  ultimate  load  carrying  capability  of  the 
face  plate  containing  the  assumed  damage  as  a  factor  of  the  initial  biaxial  loads  equal  to 


53.47. 


11  Conclusion 

A  general  predictive  methodology  for  determining  residual  strength  in  impact  damaged  com¬ 
posite  laminates  has  been  developed  and  incorporated  into  a  computer  code  designated  RE- 
STRAN  (REsidual  STRength  Analysis).  RESTRAN  can  be  used  to  analyze  composite 
structures  with  arbitrary  three-dimensional  geometry,  loading  and  support  conditions,  ma¬ 
terial  properties,  and  initial  material  and  delamination  damage.  Material  failure  modes  are 
predicted  using  a  robust  suite  of  failure  criteria  and  damage  laws.  Structural  failure  due 
to  sequential  sublaminate  buckling  of  delaminated  layers  is  also  accounted.  A  progressive 
failure  analysis  is  performed  until  ultimate  structural  failure  is  predicted  yielding  an  estimate 
of  the  residual  strength. 
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1  COMMANDER 

US  ARMY  ARDEC 
AMSTA  ARFSFT 
C  LIVECCHIA 
PICATINNY  ARSENAL  NJ 
07806-5000 

1  COMMANDER 

US  ARMY  ARDEC 
AMSTA  ASF 

PICATINNY  ARSENAL  NJ 
07806-5000 
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NO.  OF 
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1 


1 


3 


2 


1 


2 


NO.  OF 


ORGANIZATION 

COPIES 

ORGANIZATION 

COMMANDER 

1 

COMMANDER 

US  ARMY  ARDEC 

US  ARMY  ARDEC 

AMSTA  AR  QAC  T  C 

AMSTA  AR  WET 

C  PATEL 

TSACHAR 

PICATINNY  ARSENAL  NJ 

BLDG  172 

07806-5000 

COMMANDER 

PICATINNY  ARSENAL  NJ 
07806-5000 

US  ARMY  ARDEC 

9 

COMMANDER 

AMSTA  ARM 

US  ARMY  ARDEC 

D  DEMELLA 

AMSTA  ARCCHB 

PICATINNY  ARSENAL  NJ 

P  DONADIA 

07806-5000 

F  DONLON 

P  VALENTI 

COMMANDER 

C  KNUTSON 

US  ARMY  ARDEC 

G  EUSTICE 

AMSTA  ARFSA 

S  PATEL 

AWARNASH 

G  WAGNECZ 

BMACHAK 

RSAYER 

M  CHIEFA 

F CHANG 

PICATINNY  ARSENAL  NJ 

PICATINNY  ARSENAL  NJ 

07806-5000 

07806-5000 

COMMANDER 

6 

COMMANDER 

US  ARMY  ARDEC 

US  ARMY  ARDEC 

AMSTA  ARFSPG 

AMSTA  AR  CCL 

M  SCHIKSNIS 

F  PUZYCKI 

D  CARLUCCI 

R  MCHUGH 

PICATINNY  ARSENAL  NJ 

D  CONWAY 

07806-5000 

EJAROSZEWSKI 

R  SCHLENNER 

COMMANDER 

MCLUNE 

US  ARMY  ARDEC 

PICATINNY  ARSENAL  NJ 

AMSTAARFSPA 

P  KISATSKY 

07806-5000 

PICATINNY  ARSENAL  NJ 

5 

PM  SAD  ARM 

07806-5000 

SFAE  GCSS  SD 

COLB  ELLIS 

COMMANDER 

M  DEVINE 

US  ARMY  ARDEC 

W  DEMASSI 

AMSTA  ARCCHC 

J  PRITCHARD 

HCHANIN 

SHROWNAK 

S  CHICO 

PICATINNY  ARSENAL  NJ 

PICATINNY  ARSENAL  NJ 
07806-5000 

07806-5000 

1 

US  ARMY  ARDEC 

COMMANDER 

INTELLIGENCE  SPECIALIST 

US  ARMY  ARDEC 

AMSTA  ARWELF 

AMSTA  AR  QAC  T 

M  GUERRIERE 

D  RIGOGLIOSO 

PICATINNY  ARSENAL  NJ 

PICATINNY  ARSENAL  NJ 
07806-5000 

07806-5000 
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12 


1 


1 


1 


ORGANIZATION 


NO.  OF 
COPIES 


PEO  FIELD  ARTILLERY  SYS  3 

SFAE  FAS  PM 

H  GOLDMAN 

T  MCWILLIAMS 

PICATINNY  ARSENAL  NJ 

07806-5000 

PM  TMAS 
SFAE  GSSC  TMA 

R  MORRIS  1 

CKIMKER 

D  GUZIEWICZ 

EKOPACZ 

R  ROESER 

R DARCY 

R  KOWALSKI 

RMCDANOLDS  1 

L  D  ULISSE 
C  ROLLER 
J  MCGREEN 
B  PATTER 

PICATINNY  ARSENAL  NJ 
07806-5000 

1 

COMMANDER 
US  ARMYARDEC 
AMSTA  AR  WEA 
J  BRESCIA 

PICATINNY  ARSENAL  NJ 
07806-5000 

COMMANDER  2 

US  ARMY  ARDEC 

PRODUCTION  BASE 

MODERN  ACTY 

AMSMC  PBMK 

PICATINNY  ARSENAL  NJ 

07806-5000 

COMMANDER 
US  ARMY  TACOM 

PM  ABRAMS  1 

SFAE  ASM  AB 

6501  ELEVEN  MILE  RD 

WARREN  MI  48397-5000 

COMMANDER 
US  ARMY  TACOM 
AMSTA  SF 

WARREN  MI  48397-5000 


ORGANIZATION 


COMMANDER 

US  ARMY  TACOM 

PM  TACTICAL  VEHICLES 

SFAE  TVL 

SFAE  TVM 

SFAE  TVH 

6501  ELEVEN  MILE  RD 
WARREN  MI  48397-5000 

COMMANDER 
US  ARMY  TACOM 
PM  BFVS 
SFAE  ASM  BV 
6501  ELEVEN  MILE  RD 
WARREN  MI  48397-5000 

COMMANDER 

US  ARMY  TACOM 

PMAFAS 

SFAE  ASM  AF 

6501  ELEVEN  MILE  RD 

WARREN  MI  48397-5000 

COMMANDER 
US  ARMY  TACOM 
PM  RDT&E 
SFAEGCSS  WAB 
J  GODELL 

6501  ELEVEN  MILE  RD 
WARREN  MI  48397-5000 

COMMANDER 
US  ARMY  TACOM 
PM  SURV  SYS 
SFAE  ASM  SS 
T  DEAN 

SFAEGCSSWGSIM 
D COCHRAN 
6501  ELEVEN  MILE  RD 
WARREN  MI  48397-5000 

US  ARMY  CERL 
RLAMPO 

2902  NEWMARK  DR 
CHAMPAIGN  IL  61822 
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COPIES 
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ORGANIZATION 

1 

COMMANDER 

US  ARMYTACOM 

PM  SURVTVABLE  SYSTEMS 

SFAE  GCSS  W  GSI H 

MRYZYI 

6501  ELEVEN  MILE  RD 

WARREN  MI  48397-5000 

14 

COMMANDER 

US  ARMY  TACOM 
AMSTA  TRR 

R  MCCLELLAND 

D  THOMAS 

J BENNETT 

D  HANSEN 

AMSTA  JSK 

1 

COMMANDER 

US  ARMYTACOM 

PMBFV 

SFAE  GCSS  W  BV 

S  DAVIS 

6501  ELEVEN  MILE  RD 

WARREN  MI  48397-5000 

S  GOODMAN 

J FLORENCE 

KIYER 

D  TEMPLETON 

A  SCHUMACHER 
AMSTA  TRD 

D  OSTBERG 

L  HINOIOSA 

1 

COMMANDER 

US  ARMYTACOM 

CHIEF  ABRAMS  TESTING 

SFAE  GCSS  W  AB  QT 
TKRASKIEWICZ 

6501  ELEVEN  MILE  RD 

BRAJU 

AMSTA  CS  SF 

H  HUTCHINSON 

F  SCHWARZ 

WARREN  MI  48397-5000 

WARREN  MI  48397-5000  \A 

1  COMMANDER 
WATERVLIET  ARSENAL 
SMCWV  QAE  Q 

B  VANINA 
BLDG  44 

WATERVLIET  NY  12 1 89-4050 

2  TSM  ABRAMS 
ATZKTS 

S  JABURG 
W  MEINSHAUSEN 
FT  KNOX  KY  40121 

3  ARMOR  SCHOOL 
ATZKTD 
RBAUEN 
JBERG 

A  POMEY  2 

FT  KNOX  KY  40121 


BENET  LABORATORIES 
AMSTA  AR  CCB 
R  FISCELLA 
MSOJA 
EKATHE 
M  SCAVULO 
G  SPENCER 
P  WHEELER 
SKRUPSKI 
J  VASILAKIS 
G  FRIAR 
R  HASENBEIN 
AMSTA  CCB  R 
S  SOPOK 
E HYLAND 
D  CRAYON 
R  DILLON 

WATERVLIET  NY  12189-4050 

HQ  IOC  TANK 
AMMUNITION  TEAM 
AMSIO  SMT 
R  CRAWFORD 
W  HARRIS 

ROCK  ISLAND  IL  61299-6000 


NO.  OF 

COPIES  ORGANIZATION 


NO.  OF 

COPIES  ORGANIZATION 


2  COMMANDER 

USARMYAMCOM 
AVIATION  APPLIED  TECH  DIR 
J  SCHUCK 

FT  EUSTIS  VA  23604-5577 

1  DIRECTOR 
USARMYAMCOM 
SFAE  AV  RAM  TV 
D  CALDWELL 
BLDG  5300 

REDSTONE  ARSENAL  AL 
35898 

2  US  ARMY  CORPS  OF  ENGINEERS 
CERD  C 

T  LIU 
CEW  ET 
T  TAN 

20  MASS  AVE  NW 
WASHINGTON  DC  203 14 

1  US  ARMY  COLD  REGIONS 

RSCH  &  ENGRNG  LAB 
PDUTTA 
72  LYME  RD 
HANOVER  NH  03755 

1  SYSTEM  MANAGER  ABRAMS 

ATZKTS 
LTCJHNUNN 
BLDG  1002  RM  110 
FT  KNOX  KY  40121 

1  USA  SBCCOM  PM  SOLDIER  SPT 
AMSSB  PM  RSS  A 

J  CONNORS 
KANSAS  ST 
NATICK  MA  01760-5057 

2  USA  SBCCOM 
MATERIAL  SCIENCE  TEAM 
AMSSB  RSS 

J  HERBERT 
M  SENNETT 
KANSAS  ST 
NATICK  MA  01760-5057 


2  OFC  OF  NAVAL  RESEARCH 
D  SIEGEL  CODE  351 
J  KELLY 

800  N  QUINCY  ST 
ARLINGTON  VA  22217-5660 

1  NAVAL  SURFACE  WARFARE  CTR 

DAHLGREN  DIV  CODE  G06 
DAHLGREN  VA  22448 

1  NAVAL  SURFACE  WARFARE  CTR 

TECH  LIBRARY  CODE  323 
17320  DAHLGREN  RD 
DAHLGREN  VA  22448 

1  NAVAL  SURFACE  WARFARE  CTR 

CRANE  DIVISION 
M  JOHNSON  CODE  20H4 
LOUISVILLE  KY  40214-5245 

8  DIRECTOR 

US  ARMY  NATIONAL  GROUND 
INTELLIGENCE  CTR 
D  LEITER 
MHOLTUS 
M  WOLFE 
S  MINGLEDORF 
J  GASTON 

W  GSTATTENBAUER 
R  WARNER 
J  CRIDER 

220  SEVENTH  ST  NE 
CHARLOTTESVILLE  VA  22091 

2  NAVAL  SURFACE  WARFARE  CTR 

U  SORATHIA 
C  WILLIAMS  CD  6551 
9500  MACARTHUR  BLVD 
WEST  BETHESDA  MD  20817 

2  COMMANDER 

NAVAL  SURFACE  WARFARE  CTR 
CARDEROCK  DIVISION 
R  PETERSON  CODE  2020 
M  CRITCHFIELD  CODE  1730 
BETHESDA  MD  20084 
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8  US  ARMY  SBCCOM 
SOLDIER  SYSTEMS  CENTER 
BALLISTICS  TEAM 

J  WARD 
WZUKAS 
P  CUNNIFF 
J  SONG 

MARINE  CORPS  TEAM 

J  MACKIEWICZ 

BUS  AREA  ADVOCACY  TEAM 

W  HASKELL 

AMSSB  RCP  SS 

W  NYKVIST 

S  BEAUDOIN 

KANSAS  ST 

NATICK  MA  01760-5019 

9  US  ARMY  RESEARCH  OFC 
A  CROWSON 

J CHANDRA 
H EVERETT 
J  PRATER 
R  SINGLETON 
G  ANDERSON 
D  STEPP 
D  KISEROW 
J CHANG 
PO  BOX  12211 

RESEARCH  TRIANGLE  PARK  NC 
27709-2211 

8  NAVAL  SURFACE  WARFARE  CTR 
J  FRANCIS  CODE  G30 
D  WILSON  CODE  G32 
R  D  COOPER  CODE  G32 
J  FRAYSSE  CODE  G33 
E  ROWE  CODE  G33 
T  DURAN  CODE  G33 
L  DE  SIMONE  CODE  G33 
R  HUBBARD  CODE  G33 
DAHLGREN  VA  22448 

1  NAVAL  SEA  SYSTEMS  CMD 

D  LIESE 

2531  JEFFERSON  DAVIS  HWY 
ARLINGTON  VA  22242-5 1 60 

1  NAVAL  SURFACE  WARFARE  CTR 

M  LACY  CODE  B02 
17320  DAHLGREN  RD 
DAHLGREN  VA  22448 


2  NAVAL  SURFACE  WARFARE  CTR 

CARDEROCK  DIVISION 
R  CRANE  CODE  2802 
C  WILLIAMS  CODE  6553 
3A  LEGGETT  CIR 
BETHESDA  MD  20054-5000 

1  EXPEDITIONARY  WARFARE 

DIVN85 
FSHOUP 

2000  NAVY  PENTAGON 
WASHINGTON  DC  20350-2000 

1  AFRL  MLBC 

2941  P  ST  RM  136 

WRIGHT  PATTERSON  AFB  OH 

45433-7750 

1  AFRLMLSS 
R THOMSON 

2179  12THSTRM  122 
WRIGHT  PATTERSON  AFB  OH 
45433-7718 

2  AFRL 

F  ABRAMS 

J  BROWN 

BLDG  653 

2977  P  ST  STE  6 

WRIGHT  PATTERSON  AFB  OH 

45433-7739 

1  WATERWAYS  EXPERIMENT 

D  SCOTT 

3909  HALLS  FERRY  RD  SC  C 
VICKSBURG  MS  39180 

5  DIRECTOR 

LLNL 

R  CHRISTENSEN 
S  DETERESA 
F  MAGNESS 
M  FINGER  MS  313 
M  MURPHY  L  282 
PO  BOX  808 
LIVERMORE  CA  94550 

1  AFRL  MLS  OL 

L  COULTER 
7278  4TH  ST 
BLDG  100  BAY  D 
HILL  AFB  UT  84056-5205 
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COPIES  ORGANIZATION 


NO.  OF 

COPIES  ORGANIZATION 


1  OSD 

JOINT  CCD  TEST  FORCE 
OSD  JCCD 
R  WILLIAMS 
3909  HALLS  FERRY  RD 
VICKSBURG  MS  29180-6199 

1  DEFENSE  NUCLEAR  AGENCY 
INNOVATIVE  CONCEPTS  DIV 
6801  TELEGRAPH  RD 
ALEXANDRIA  VA  22310-3398 

3  DARPA 

M  VANFOSSEN 
SWAX 

L  CHRISTODOULOU 
3701  N  FAIRFAX  DR 
ARLINGTON  VA  22203-1714 

2  SERDP  PROGRAM  OFC 
PMP2 

C  PELLERIN 
B  SMITH 

901  N  STUART  ST  STE  303 
ARLINGTON  VA  22203 

1  FAA 

MILHDBK  17  CHAIR 
LILCEWICZ 
1601  LIND  AVE  SW 
ANM  115N 
RESTONVA  98055 

1  US  DEPT  OF  ENERGY 

OFC  OF  ENVIRONMENTAL 
MANAGEMENT 
P  RITZCOVAN 
19901  GERMANTOWN  RD 
GERMANTOWN  MD  20874-1928 

1  DIRECTOR 

LLNL 

F  ADDESSIO  MSB216 

PO  BOX  1633 

LOS  ALAMOS  NM  87545 

1  OAK  RIDGE  NATIONAL 

LABORATORY 
RM  DAVIS 
PO  BOX  2008 

OAK  RIDGE  TN  37831-6195 


3  DIRECTOR 

SANDIA  NATIONAL  LABS 

APPLIED  MECHANICS  DEPT 

MS  9042 

JHANDROCK 

YRKAN 

J  LAUFFER 

PO  BOX  969 

LIVERMORE  CA  94551-0969 

1  OAK  RIDGE  NATIONAL 

LABORATORY 
CEBERLEMS  8048 
PO  BOX  2008 
OAK  RIDGE  TN  37831 

1  OAK  RIDGE  NATIONAL 

LABORATORY 
CD  WARREN  MS  8039 
PO  BOX  2008 
OAK  RIDGE  TN  37831 

5  NIST 

JDUNKERS 

M  VANLANDINGHAM  MS  8621 
J  CHIN  MS  8621 
J  MARTIN  MS  8621 
D  DUTHINH  MS  8611 
100  BUREAU  DR 
GAITHERSBURG  MD  20899 

1  HYDROGEOLOGIC  INC 

SERDP  ESTCP  SPT  OFC 
S  WALSH 

1 155  HERNDON  PKWY  STE  900 
HERNDON  VA  20170 

3  NASA  LANGLEY  RSCH  CTR 

AMSRL  VS 
WELBERMS  266 
F  BARTLETT  JR  MS  266 
G  FARLEY  MS  266 
HAMPTON  VA  23681-0001 

1  NASA  LANGLEY  RSCH  CTR 

T  GATES  MS  188E 
HAMPTON  VA  23661-3400 

1  FHWA 

EMUNLEY 

6300  GEORGETOWN  PIKE 
MCLEAN  VA  22101 
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3  CYTEC  FIBERITE 

R DUNNE 
DKOHLI 
RMAYHEW 
1300  REVOLUTION  ST 
HAVRE  DE  GRACE  MD  21078 

1  USDOT  FEDERAL  RAILRD 

M  FATEH  RDV  31 
WASHINGTON  DC  20590 

1  MARINE  CORPS 

INTLLGNC  ACTVTY 
D  KOSITZKE 

3300  RUSSELL  RD  STE  250 
QUANTICO  VA  22134-501 1 

1  DIRECTOR 

NATIONAL  GRND  INTLLGNC  CTR 

IANGTMT 

220  SEVENTH  ST  NE 

CHARLOTTESVILLE  VA 

22902-5396 

1  SIOUX  MFG 
BKRIEL 
PO  BOX  400 

FT  TOTTEN  ND  58335 

2  3TEX  CORPORATION 
A  BOGDANOVICH 

J  SINGLETARY 
109  MACKENAN  DR 
CARY  NC  275 11 

1  3M  CORPORATION 

J  SKILDUM 

3M  CENTER  BLDG  60  IN  01 
ST  PAUL  MN  55144-1000 

1  DIRECTOR 

DEFENSE  INTLLGNC  AGNCY 
TA  5 

K  CRELLING 
WASHINGTON  DC  20310 

1  ADVANCED  GLASS  FIBER  YARNS 

T  COLLINS 

281  SPRING  RUN  LANE  STE  A 
DOWNINGTON  PA  19335 


1  COMPOSITE  MATERIALS  INC 
D  SHORTT 
19105  63  AVENE 
PO  BOX  25 

ARLINGTON  WA  98223 

1  JPS  GLASS 

L  CARTER 
PO  BOX  260 
SLATER  RD 
SLATER  SC  29683 

1  COMPOSITE  MATERIALS  INC 

R  HOLLAND 
11  JEWEL  CT 
ORINDA  CA  94563 

1  COMPOSITE  MATERIALS  INC 
C  RILEY 

14530  S  ANSON  AVE 
SANTA  FE  SPRINGS  CA  90670 

2  SIMULA 

J  COLTMAN 
RHUYETT 
10016  S51ST  ST 
PHOENIX  AZ  85044 

2  PROTECTION  MATERIALS  INC 

M  MILLER 
F  CRILLEY 
14000  NW  58  CT 
MIAMI  LAKES  FL  33014 

2  FOSTER  MILLER 

M  ROYLANCE 
WZUKAS 
195  BEAR  HILL  RD 
WALTHAM  MA  02354-1 196 

1  ROM  DEVELOPMENT  CORP 

ROMEARA 
136  SWINEBURNE  ROW 
BRICK  MARKET  PLACE 
NEWPORT  RI 02840 
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2  TEXTRON  SYSTEMS 

T  FOLTZ 
M TREASURE 
1449  MIDDLESEX  ST 
LOWELL  MA  01851 

1  O  GARA  HESS  &  EISENHARDT 
M  GILLESPIE 

9113  LESAINT  DR 
FAIRFIELD  OH  450 14 

2  MILLIKEN  RSCH  CORP 
H  KUHN 

M  MACLEOD 
PO  BOX  1926 
SPARTANBURG  SC  29303 

1  CONNEAUGHT  INDUSTRIES  INC 

J  SANTOS 
PO  BOX  1425 
COVENTRY  RI 02816 

1  BATTELLE  NATICK  OPNS 

BHALPIN 

209  W  CENTRAL  ST  STE  302 
NATICK  MA  01760 

1  ARMTEC  DEFENSE  PRODUCTS 

S  DYER 
85  901  AVE  53 
PO  BOX  848 
COACHELLA  CA  92236 

1  NATIONAL  COMPOSITE  CENTER 
T  CORDELL 

2000  COMPOSITE  DR 
KETTERING  OH  45420 

3  PACIFIC  NORTHWEST  LAB 
M  SMITH 

G  VAN  ARSDALE 
R  SHIPPELL 
PO  BOX  999 
RICHLAND  WA  99352 

2  AMOCO  PERFORMANCE 
PRODUCTS 

M  MICHNO  JR 
J  BANISAUKAS 
4500  MCGINNIS  FERRY  RD 
ALPHARETTA  GA  30202-3944 


8  ALLIANT  TECHSYSTEMS  INC 
C  CANDLAND  MN1 1  2830 
C  AAKHUS  MN1 1  2830 
B  SEE  MN11  2439 
N  VLAHAKUS  MN1 1  2145 
RDOHRNMN11  2830 
S  HAGLUND  MN1 1  2439 
M  HISSONG  MN1 1  2830 
D  KAMDAR  MN1 1  2830 
600  SECOND  ST  NE 
HOPKINS  MN  55343-8367 

1  SAIC 

M  PALMER 

1410  SPRING  HILL  RD  STE  400 
MS  SH4 5 

MCLEAN  VA  22 102 

1  SAIC 

G  CHRYSSOMALLIS 
3800  W  80TH  ST  STE  1090 
BLOOMINGTON  MN  55431 

1  AAI  CORPORATION 

T  G  STASTNY 
PO  BOX  126 

HUNT  VALLEY  MD  21030-0126 

1  APPLIED  COMPOSITES 

W  GRISCH 

333  NORTH  SIXTH  ST 
ST  CHARLES  IL  60174 

1  CUSTOM  ANALYTICAL 

ENG  SYS  INC 
A  ALEXANDER 
13000  TENSOR  LANE  NE 
FLINTSTONE  MD  21530 

3  ALLIANT  TECHSYSTEMS  INC 

J  CONDON 
ELYNAM 
J GERHARD 
WV01  16  STATE  RT  956 
PO  BOX  210 

ROCKET  CENTER  WV  26726-0210 

1  OFC  DEPUTY  UNDER  SEC  DEFNS 

J  THOMPSON 

1745  JEFFERSON  DAVIS  HWY 
CRYSTAL  SQ  4  STE  501 
ARLINGTON  VA  22202 
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1  PROJECTILE  TECHNOLOGY  INC 

515  GILES  ST 

HAVRE  DE  GRACE  MD  21078 

5  AEROJET  GEN  CORP 

D  PILLASCH 
T  COULTER 
C  FLYNN 
DRUBAREZUL 
M  GREINER 

1100  WEST  HOLLYVALE  ST 
AZUSA  CA  91702-0296 

3  HEXCEL  INC 

RBOE 

PO  BOX  18748 

SALT  LAKE  CITY  UT  84 1 1 8 

1  HERCULES  INC 

HERCULES  PLAZA 
WILMINGTON  DE  19894 

1  BRIGS  COMPANY 

J  BACKOFEN 
2668  PETERBOROUGH  ST 
HERNDON  VA  22071-2443 

1  ZERNOW  TECHNICAL  SERVICES 

L  ZERNOW 

425  W  BONITA  AVE  STE  208 
SAN  DIMAS  CA  91773 

1  GENERAL  DYNAMICS  OTS 

L  WHITMORE 
10101  NINTH  ST  NORTH 
ST  PETERSBURG  FL  33702 

3  GENERAL  DYNAMICS  OTS 

FLINCHBAUGH  DIV 
E  STEINER 
B  STEWART 
T LYNCH 
PO  BOX  127 
RED  LION  PA  17356 

1  GKN  AEROSPACE 

DOLDS 

15  STERLING  DR 
WALLINGFORD  CT  06492 


5  SIKORSKY  AIRCRAFT 
G  JACARUSO 
T  CARSTENSAN 
B  KAY 

S  GARBO  MS  S330A 
J  ADELMANN 
6900  MAIN  ST 
PO  BOX  9729 

STRATFORD  CT  06497-9729 

1  PRATT  &  WHITNEY 

C  WATSON 

400  MAIN  ST  MS  114  37 
EAST  HARTFORD  CT  06108 

1  AEROSPACE  CORP 
G  HAWKINS  M4  945 

2350  E  EL  SEGUNDO  BLVD 
EL  SEGUNDO  CA  90245 

2  CYTEC  FIBERITE 
M  LIN 

W  WEB 

1440  N  KRAEMER  BLVD 
ANAHEIM  CA  92806 

1  UDLP 

G  THOMAS 

PO  BOX  58123 

SANTA  CLARA  CA  95052 

2  UDLP 

R  BARRETT  MAIL  DROP  M53 
V  HORVATICH  MAIL  DROP  M53 
328  WBROKAWRD 
SANTA  CLARA  CA  95052-0359 

3  UDLP 

GROUND  SYSTEMS  DIVISION 
M  PEDRAZZI  MAIL  DROP  N09 
A  LEE  MAIL  DROP  Nil 
M  MACLEAN  MAIL  DROP  N06 
1205  COLEMAN  AVE 
SANTA  CLARA  CA  95052 

4  UDLP 

R  BRYNSVOLD 
P  JANKE  MS  170 
4800  EAST  RIVER  RD 
MINNEAPOLIS  MN  55421-1498 
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